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Series Preface 



Signal processing applications are now widespread. Relatively cheap consumer 
products through to the more expensive military and industrial systems extensively 
exploit this technology. This spread was initiated in the 1960s by the introduction of 
cheap digital technology to implement signal processing algorithms in real-time for 
some applications. Since that time semiconductor technology has developed rapidly 
to support the spread. In parallel, an ever increasing body of mathematical theory 
is being used to develop signal processing algorithms. The basic mathematical 
foundations, however, have been known and well understood for some time. 

Signal Processing and its Applications addresses the entire breadth and depth 
of the subject with texts that cover the theory, technology and applications of signal 
processing in its widest sense. This is reflected in the composition of the Editorial 
Board, who have interests in: 

(i) Theory - The physics of the application and the mathematics to model the 
system; 

(ii) Implementation - VLSI/ ASIC design, computer architecture, numerical 
methods, systems design methodology, and CAE; 

(iii) Applications - Speech, sonar, radar, seismic, medical, communications (both 
audio and video), guidance, navigation, remote sensing, imaging, survey, 
archiving, non-destructive and non-intrusive testing, and personal entertain- 
ment. 

Signal Processing and its Applications will typically be of most interest to post- 
graduate students, academics, and practising engineers who work in the held and 
develop signal processing applications. Some texts may also be of interest to final 
year undergraduates. 

Richard C. Green 
The Engineering Practice , 
Farnborough, UK 
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Preface 



It is hoped that this book will serve both as a text in time-series analysis and signal 
processing and as a reference book for research workers and practitioners. Time- 
series analysis and signal processing are two subjects which ought to be treated 
as one; and they are the concern of a wide range of applied disciplines includ- 
ing statistics, electrical engineering, mechanical engineering, physics, medicine and 
economics. 

The book is primarily a didactic text and, as such, it has three main aspects. 
The first aspect of the exposition is the mathematical theory which is the foundation 
of the two subjects. The book does not skimp this. The exposition begins in 
Chapters 2 and 3 with polynomial algebra and complex analysis, and it reaches 
into the middle of the book where a lengthy chapter on Fourier analysis is to be 
found. 

The second aspect of the exposition is an extensive treatment of the numerical 
analysis which is specifically related to the subjects of time-series analysis and 
signal processing but which is, usually, of a much wider applicability. This be- 
gins in earnest with the account of polynomial computation, in Chapter 4, and 
of matrix computation, in Chapter 7, and it continues unabated throughout the 
text. The computer code, which is the product of the analysis, is distributed 
evenly throughout the book, but it is also hierarchically ordered in the sense that 
computer procedures which come later often invoke their predecessors. 

The third and most important didactic aspect of the text is the exposition of 
the subjects of time-series analysis and signal processing themselves. This begins 
as soon as, in logic, it can. However, the fact that the treatment of the substantive 
aspects of the subject is delayed until the mathematical foundations are in place 
should not prevent the reader from embarking immediately upon such topics as the 
statistical analysis of time series or the theory of linear filtering. The book has been 
assembled in the expectation that it will be read backwards as well as forwards, as 
is usual with such texts. Therefore it contains extensive cross-referencing. 

The book is also intended as an accessible work of reference. The computer 
code which implements the algorithms is woven into the text so that it binds closely 
with the mathematical exposition; and this should allow the detailed workings of 
the algorithms to be understood quickly. However, the function of each of the Pascal 
procedures and the means of invoking them are described in a reference section, 
and the code of the procedures is available in electronic form on a computer disc. 

The associated disc contains the Pascal code precisely as it is printed in the 
text. An alternative code in the C language is also provided. Each procedure is 
coupled with a so-called driver, which is a small program which shows the procedure 
in action. The essential object of the driver is to demonstrate the workings of 
the procedure; but usually it fulfils the additional purpose of demonstrating some 
aspect the theory which has been set forth in the chapter in which the code of the 
procedure it to be found. It is hoped that, by using the algorithms provided in this 
book, scientists and engineers will be able to piece together reliable software tools 
tailored to their own specific needs. 
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The compact disc also contains a collection of reference material which includes 
the libraries of the computer routines and various versions of the bibliography of 
the book. The numbers in brackets which accompany the bibliographic citations 
refer to their order in the composite bibliography which is to be found on the disc. 
On the disc, there is also a bibliography which is classified by subject area. 

A preface is the appropriate place to describe the philosophy and the motiva- 
tion of the author in so far as they affect the book. A characteristic of this book, 
which may require some justification, is its heavy emphasis on the mathematical 
foundations of its subjects. There are some who regard mathematics as a burden 
which should be eased or lightened whenever possible. The opinion which is re- 
flected in the book is that a firm mathematical framework is needed in order to bear 
the weight of the practical subjects which are its principal concern. For example, 
it seems that, unless the reader is adequately appraised of the notions underlying 
Fourier analysis, then the perplexities and confusions which will inevitably arise 
will limit their ability to commit much of the theory of linear filters to memory. 
Practical mathematical results which are well-supported by theory are far more 
accessible than those which are to be found beneath piles of technological detritus. 

Another characteristic of the book which reflects a methodological opinion is 
the manner in which the computer code is presented. There are some who regard 
computer procedures merely as technological artefacts to be encapsulated in boxes 
whose contents are best left undisturbed for fear of disarranging them. An opposite 
opinion is reflected in this book. The computer code presented here should be read 
and picked to pieces before being reassembled in whichever way pleases the reader. 
In short, the computer procedures should be approached in a spirit of constructive 
play. An individual who takes such an approach in general will not be balked by 
the non-availability of a crucial procedure or by the incapacity of some large-scale 
computer program upon which they have come to rely. They will be prepared 
to make for themselves whatever tools they happen to need for their immediate 
purposes. 

The advent of the microcomputer has enabled the approach of individualist 
self-help advocated above to become a practical one. At the same time, it has 
stimulated the production of a great variety of highly competent scientific software 
which is supplied commercially. It often seems like wasted effort to do for oneself 
what can sometimes be done better by purpose-built commercial programs. Clearly, 
there are opposing forces at work here — and the issue is the perennial one of whether 
we are to be the masters or the slaves of our technology. The conflict will never be 
resolved; but a balance can be struck. This book, which aims to help the reader to 
master one of the central technologies of the latter half of this century, places most 
of its weight on one side of the scales. 

D.S.G. POLLOCK 
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CHAPTER 1 



The Methods of 
Time-Series Analysis 



The methods to be presented in this book are designed for the purpose of analysing 
series of statistical observations taken at regular intervals in time. The methods 
have a wide range of applications. We can cite astronomy [539], meteorology [444], 
seismology [491], oceanography [232], [251], communications engineering and signal 
processing [425], the control of continuous process plants [479], neurology and elec- 
troencephalography [151], [540], and economics [233]; and this list is by no means 
complete. 

The Frequency Domain and the Time Domain 

The methods apply, in the main, to what are described as stationary or non- 
evolutionary time series. Such series manifest statistical properties which are in- 
variant throughout time, so that the behaviour during one epoch is the same as it 
would be during any other. 

When we speak of a weakly stationary or covariance-stationary process, we 
have in mind a sequence of random variables y{t) = {yt',t = 0, ±1,±2, . . .}, rep- 
resenting the potential observations of the process, which have a common finite 
expected value E(y t ) = y and a set of autocovariances C(yt,y s ) = E{(y t — y)(y s — 
y)} = 7|t- s | which depend only on the temporal separation r = \t — s| of the dates 
t and s and not on their absolute values. Usually, we require of such a process 
that lim(r — > oo)^ T = 0, which is to say that the correlation between increasingly 
remote elements of the sequence tends to zero. This is a way of expressing the 
notion that the events of the past have a diminishing effect upon the present as 
they recede in time. In an appendix to the chapter, we review the definitions of 
mathematical expectations and covariances. 

There are two distinct yet broadly equivalent modes of time-series analysis 
which may be pursued. On the one hand are the time-domain methods which 
have their origin in the classical theory of correlation. Such methods deal pre- 
ponderantly with the autocovariance functions and the cross-covariance functions 
of the series, and they lead inevitably towards the construction of structural or 
parametric models of the autoregressive moving-average type for single series and 
of the transfer-function type for two or more causally related series. Many of the 
methods which are used to estimate the parameters of these models can be viewed 
as sophisticated variants of the method of linear regression. 

On the other hand are the frequency-domain methods of spectral analysis. 
These are based on an extension of the methods of Fourier analysis which originate 
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in the idea that, over a finite interval, any analytic function can be approximated, 
to whatever degree of accuracy is desired, by taking a weighted sum of sine and 
cosine functions of harmonically increasing frequencies. 

Harmonic Analysis 

The astronomers are usually given credit for being the first to apply the meth- 
ods of Fourier analysis to time series. Their endeavours could be described as the 
search for hidden periodicities within astronomical data. Typical examples were 
the attempts to uncover periodicities within the activities recorded by the Wolfer 
sunspot index — see Izenman [266] — and in the indices of luminosity of variable 
stars. 

The relevant methods were developed over a long period of time. Lagrange 
[306] suggested methods for detecting hidden periodicities in 1772 and 1778. The 
Dutchman Buijs-Ballot [86] propounded effective computational procedures for the 
statistical analysis of astronomical data in 1847. However, we should probably 
credit Sir Arthur Schuster [444], who in 1889 propounded the technique of periodo- 
gram analysis, with being the progenitor of the modern methods for analysing time 
series in the frequency domain. 

In essence, these frequency-domain methods envisaged a model underlying the 
observations which takes the form of 

y(t) = Y Pi cosiujt - 0j) + e{t) 

(L1) A r 

= 2 ^ \ a i cos (ujt) + Pj sm (ujt)} + e(t), 
i 



where ay = p 3 cos 0 3 and P 3 = p 3 sin (i 3 , and where e(t) is a sequence of indepen- 
dently and identically distributed random variables which we call a white-noise 
process. Thus the model depicts the series y(t) as a weighted sum of perfectly 
regular periodic components upon which is superimposed a random component. 

The factor pj = yj (a? + /3|) is called the amplitude of the jth periodic com- 
ponent, and it indicates the importance of that component within the sum. Since 
the variance of a cosine function, which is also called its mean-square deviation, is 
just one half, and since cosine functions at different frequencies are uncorrelated, 
it follows that the variance of y(t) is expressible as V{y(t)} = \ P 2 j + a e where 
<7g = 14{e(f)} is the variance of the noise. 

The periodogram is simply a device for determining how much of the variance 
of y(t) is attributable to any given harmonic component. Its value at c o 3 = 2nj/T, 
calculated from a sample yo, - ■ ■ ,Vt- i comprising T observations on y(t), is given 



by 


, 2 




I (ui)=T 


(1.2) 


T 



Y Vt cos (wjt) | + | Y yt sin ( w ^) | 
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Figure 1.1. The graph of a sine function. 




Figure 1.2. Graph of a sine function with small random fluctuations superimposed. 



If y(t) does indeed comprise only a finite number of well-defined harmonic compo- 
nents, then it can be shown that 2 1(ujj)/T is a consistent estimator of p^j in the 
sense that it converges to the latter in probability as the size T of the sample of 
the observations on y{t ) increases. 

The process by which the ordinates of the periodogram converge upon the 
squared values of the harmonic amplitudes was well expressed by Yule [539] in a 
seminal article of 1927: 

If we take a curve representing a simple harmonic function of time, and 
superpose on the ordinates small random errors, the only effect is to make 
the graph somewhat irregular, leaving the suggestion of periodicity still 
clear to the eye (see Figures 1.1 and 1.2). If the errors are increased in 
magnitude, the graph becomes more irregular, the suggestion of periodic- 
ity more obscure, and we have only sufficiently to increase the errors to 
mask completely any appearance of periodicity. But, however large the 
errors, periodogram analysis is applicable to such a curve, and, given a 
sufficient number of periods, should yield a close approximation to the 
period and amplitude of the underlying harmonic function. 
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Figure 1.3. Wolfer’s sunspot numbers 1749-1924. 



We should not quote this passage without mentioning that Yule proceeded to 
question whether the hypothesis underlying periodogram analysis, which postulates 
the equation under (1.1), was an appropriate hypothesis for all cases. 

A highly successful application of periodogram analysis was that of Whittaker 
and Robinson [515] who, in 1924, showed that the series recording the brightness or 
magnitude of the star T. Ursa Major over 600 days could be fitted almost exactly by 
the sum of two harmonic functions with periods of 24 and 29 days. This led to the 
suggestion that what was being observed was actually a two-star system wherein the 
larger star periodically masked the smaller, brighter star. Somewhat less successful 
were the attempts of Arthur Schuster himself [445] in 1906 to substantiate the claim 
that there is an 11-year cycle in the activity recorded by the Wolfer sunspot index 
(see Figure 1.3). 

Other applications of the method of periodogram analysis were even less suc- 
cessful; and one application which was a significant failure was its use by William 
Beveridge [51], [52] in 1921 and 1922 to analyse a long series of European wheat 
prices. The periodogram of this data had so many peaks that at least twenty 
possible hidden periodicities could be picked out, and this seemed to be many 
more than could be accounted for by plausible explanations within the realm of 
economic history. Such experiences seemed to point to the inappropriateness to 
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economic circumstances of a model containing perfectly regular cycles. A classic 
expression of disbelief was made by Slutsky [468] in another article of 1927 : 

Suppose we are inclined to believe in the reality of the strict periodicity of 
the business cycle, such, for example, as the eight-year period postulated 
by Moore [352]. Then we should encounter another difficulty. Wherein 
lies the source of this regularity? What is the mechanism of causality 
which, decade after decade, reproduces the same sinusoidal wave which 
rises and falls on the surface of the social ocean with the regularity of day 
and night? 

Autoregressive and Moving- Average Models 

The next major episode in the history of the development of time-series analysis 
took place in the time domain, and it began with the two articles of 1927 by Yule 
[539] and Slutsky [468] from which we have already quoted. In both articles, we 
find a rejection of the model with deterministic harmonic components in favour 
of models more firmly rooted in the notion of random causes. In a wonderfully 
figurative exposition, Yule invited his readers to imagine a pendulum attached to 
a recording device and left to swing. Then any deviations from perfectly harmonic 
motion which might be recorded must be the result of errors of observation which 
could be all but eliminated if a long sequence of observations were subjected to a 
periodogram analysis. Next, Yule enjoined the reader to imagine that the regular 
swing of the pendulum is interrupted by small boys who get into the room and 
start pelting the pendulum with peas sometimes from one side and sometimes from 
the other. The motion is now affected not by superposed fluctuations but by true 
disturbances. 

In this example, Yule contrives a perfect analogy for the autoregressive time- 
series model. To explain the analogy, let us begin by considering a homogeneous 
second-order difference equation of the form 

(1.3) y(t) = - 1) + </> 2 2 /(f - 2). 

Given the initial values j/_i and y~ 2 , this equation can be used recursively to 
generate an ensuing sequence {yo,yi, ■ ■ ■}■ This sequence will show a regular 
pattern of behaviour whose nature depends on the parameters <f> \ and <f> 2 - If 
these parameters are such that the roots of the quadratic equation z 2 — (j)\z — 
(j >2 = 0 are complex and less than unity in modulus, then the sequence of 
values will show a damped sinusoidal behaviour just as a clock pendulum will 
which is left to swing without the assistance of the falling weights. In fact, in 
such a case, the general solution to the difference equation will take the form 
of 

(1.4) y(t) = ap* cos (cut — 9), 

where the modulus p, which has a value between 0 and 1, is now the damp- 
ing factor which is responsible for the attenuation of the swing as the time t 
elapses. 
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Figure 1.4. A series generated by Yule’s equation 
y{t) = 1.343 y(t - 1) - 0.655 y{t - 2) + e(t). 




0 10 20 30 40 50 60 70 80 90 

Figure 1.5. A series generated by the equation 
y{t) = 1.576 y(t - 1) - 0.903 y(t - 2) + e(t). 



The autoregressive model which Yule was proposing takes the form of 
(1.5) y{t) = c/)iy(t - 1) + (j) 2 y(t - 2) + e(t), 

where e(t) is, once more, a white-noise sequence. Now, instead of masking the reg- 
ular periodicity of the pendulum, the white noise has actually become the engine 
which drives the pendulum by striking it randomly in one direction and another. Its 
haphazard influence has replaced the steady force of the falling weights. Neverthe- 
less, the pendulum will still manifest a deceptively regular motion which is liable, 
if the sequence of observations is short and contains insufficient contrary evidence, 
to be misinterpreted as the effect of an underlying mechanism. 

In his article of 1927, Yule attempted to explain the Wolfer index in terms 
of the second-order autoregressive model of equation (1.5). From the empirical 
autocovariances of the sample represented in Figure 1.3, he estimated the values 
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<f > i = 1.343 and </> 2 = —0.655. The general solution of the corresponding homoge- 
neous difference equation has a damping factor of p = 0.809 and an angular velocity 
of to = 33.96 degrees. The angular velocity indicates a period of 10.6 years which 
is a little shorter than the 11-year period obtained by Schuster in his periodogram 
analysis of the same data. In Figure 1.4, we show a series which has been gen- 
erated artificially from Yule’s equation, which may be compared with a series, in 
Figure 1.5, generated by the equation y(t) = 1.576 y(t — 1) — 0.903 y(t — 2) + e(t). 
The homogeneous difference equation which corresponds to the latter has the same 
value of u> as before. Its damping factor has the value p = 0.95, and this increase 
accounts for the greater regularity of the second series. 

Neither of our two series accurately mimics the sunspot index; although the 
second series seems closer to it than the series generated by Yule’s equation. An 
obvious feature of the sunspot index which is not shared by the artificial series is the 
fact that the numbers are constrained to be nonnegative. To relieve this constraint, 
we might apply to Wolf’s numbers yt a transformation of the form log (yt + A) or 
of the more general form (y t + A) K_1 , such as has been advocated by Box and Cox 
[69]. A transformed series could be more closely mimicked. 

The contributions to time-series analysis made by Yule [539] and Slutsky [468] 
in 1927 were complementary: in fact, the two authors grasped opposite ends of 
the same pole. For ten years, Slutsky’s paper was available only in its original 
Russian version; but its contents became widely known within a much shorter 
period. 

Slutsky posed the same question as did Yule, and in much the same man- 
ner. Was it possible, he asked, that a definite structure of a connection between 
chaotically random elements could form them into a system of more or less regular 
waves? Slutsky proceeded to demonstrate this possibility by methods which were 
partly analytic and partly inductive. He discriminated between coherent series 
whose elements were serially correlated and incoherent or purely random series of 
the sort which we have described as white noise. As to the coherent series, he 
declared that 

their origin may be extremely varied, but it seems probable that an espe- 
cially prominent role is played in nature by the process of moving summa- 
tion with weights of one kind or another; by this process coherent series 
are obtained from other coherent series or from incoherent series. 

By taking, as his basis, a purely random series obtained by the People’s Com- 
missariat of Finance in drawing the numbers of a government lottery loan, and 
by repeatedly taking moving summations, Slutsky was able to generate a series 
which closely mimicked an index, of a distinctly undulatory nature, of the English 
business cycle from 1855 to 1877. 

The general form of Slutsky’s moving summation can be expressed by writing 

(1.6) y(t) = p 0 e(t) + ii^e{t - 1) H b p q e(t - q), 

where e(t) is a white-noise process. This is nowadays called a gth-order moving- 
average model, and it is readily compared to an autoregressive model of the sort 
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depicted under (1.5). The more general pth-order autoregressive model can be 
expressed by writing 



Thus, whereas the autoregressive process depends upon a linear combination of the 
function y(t) with its own lagged values, the moving-average process depends upon 
a similar combination of the function e(t) with its lagged values. The affinity of the 
two sorts of process is further confirmed when it is recognised that an autoregressive 
process of finite order is equivalent to a moving-average process of infinite order 
and that, conversely, a finite-order moving-average process is just an infinite-order 
autoregressive process. 

Generalised Harmonic Analysis 

The next step to be taken in the development of the theory of time series was 
to generalise the traditional method of periodogram analysis in such a way as to 
overcome the problems which arise when the model depicted under (1.1) is clearly 
inappropriate. 

At first sight, it would not seem possible to describe a covariance-stationary 
process, whose only regularities are statistical ones, as a linear combination of 
perfectly regular periodic components. However, any difficulties which we might 
envisage can be overcome if we are prepared to accept a description which is in 
terms of a nondenumerable infinity of periodic components. Thus, on replacing the 
so-called Fourier sum within equation (1.1) by a Fourier integral, and by deleting the 
term e(f), whose effect is now absorbed by the integrand, we obtain an expression 
in the form of 



can be no presumption that the functions A(iv) and B(ui) are continuous. As it 
stands, this expression is devoid of any statistical interpretation. Moreover, if we 
are talking of only a single realisation of the process y(t), then the generalised 
functions A{oj) and B{to) will reflect the unique peculiarities of that realisation and 
will not be amenable to any systematic description. 

However, a fruitful interpretation can be given to these functions if we consider 
the observable sequence y(t) = {y t ; t = 0, ±1, ±2, . . .} to be a particular realisation 
which has been drawn from an infinite population representing all possible reali- 
sations of the process. For, if this population is subject to statistical regularities, 
then it is reasonable to regard dA{u > ) and dB(cu) as mutually uncorrelated random 
variables with well-defined distributions which depend upon the parameters of the 
population. 

We may therefore assume that, for any value of u, 



(1.7) 



a 0 y(t) + aiy(t - 1) H b a p y(t - p) = e(t). 



(1.8) 




Here we write dA(u>) and dB(uj) rather than a(oj)dco and P(uj)doj because there 



(1.9) 



E{dA(uj)} = E{dB(u)} = 0 and 
E{dA(u)dB(u)} = 0. 
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Figure 1.6. The spectrum of the process y{t) = 1.343y(£ — 1) — 0.655j/(£ — 
2) + e(£) which generated the series in Figure 1.4. A series of a more regular 
nature would be generated if the spectrum were more narrowly concentrated 
around its modal value. 



Moreover, to express the discontinuous nature of the generalised functions, we as- 
sume that, for any two values u> and A in their domain, we have 

(1.10) E{dA(u})dA(\)} = E{dB(to)dB{\)} = 0, 

which means that A(oj) and B(u) are stochastic processes — indexed on the 
frequency parameter u rather than on time — which are uncorrelated in non- 
overlapping intervals. Finally, we assume that dA(uj) and dB(to) have a common 
variance so that 



(1.11) V{dA{u)} = V{dB(u;)} = dG{oj). 

Given the assumption of the mutual uncorrelatedness of dA{uS) and dB(u> ), it 
therefore follows from (1.8) that the variance of y(t) is expressible as 



V{y(t)}= / [cos 2 (uit)V {dA(co)} + sin 2 (ujt)V{dB(uj)}] 



( 1 . 12 ) 



= f dG(u). 
Jo 



The function G(w), which is called the spectral distribution, tells us how much of 
the variance is attributable to the periodic components whose frequencies range 
continuously from 0 to u. If none of these components contributes more than 
an infinitesimal amount to the total variance, then the function G{oj) is absolutely 
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continuous, and we can write dG(ui) = g(u>)duj under the integral of equation (1.11). 
The new function g(u)), which is called the spectral density function or the spectrum, 
is directly analogous to the function expressing the squared amplitude which is 
associated with each component in the simple harmonic model discussed in our 
earlier sections. Figure 1.6 provides an an example of a spectral density function. 

Smoothing the Periodogram 

It might be imagined that there is little hope of obtaining worthwhile estimates 
of the parameters of the population from which the single available realisation y(t) 
has been drawn. However, provided that y(t) is a stationary process, and provided 
that the statistical dependencies between widely separated elements are weak, the 
single realisation contains all the information which is necessary for the estimation 
of the spectral density function. In fact, a modified version of the traditional 
periodogram analysis is sufficient for the purpose of estimating the spectral density. 

In some respects, the problems posed by the estimation of the spectral density 
are similar to those posed by the estimation of a continuous probability density func- 
tion of unknown functional form. It is fruitless to attempt directly to estimate the 
ordinates of such a function. Instead, we might set about our task by constructing a 
histogram or bar chart to show the relative frequencies with which the observations 
that have been drawn from the distribution fall within broad intervals. Then, by 
passing a curve through the mid points of the tops of the bars, we could construct 
an envelope that might approximate to the sought-after density function. A more 
sophisticated estimation procedure would not group the observations into the fixed 
intervals of a histogram; instead it would record the number of observations falling 
within a moving interval. Moreover, a consistent method of estimation, which aims 
at converging upon the true function as the number of observations increases, would 
vary the width of the moving interval with the size of the sample, diminishing it 
sufficiently slowly as the sample size increases for the number of sample points 
falling within any interval to increase without bound. 

A common method for estimating the spectral density is very similar to the 
one which we have described for estimating a probability density function. Instead 
of being based on raw sample observations as is the method of density-function 
estimation, it is based upon the ordinates of a periodogram which has been fitted 
to the observations on y{t). This procedure for spectral estimation is therefore 
called smoothing the periodogram. 

A disadvantage of the procedure, which for many years inhibited its widespread 
use, lies in the fact that calculating the periodogram by what would seem to be 
the obvious methods be can be vastly time-consuming. Indeed, it was not until the 
mid 1960s that wholly practical computational methods were developed. 

The Equivalence of the Two Domains 

It is remarkable that such a simple technique as smoothing the periodogram 
should provide a theoretical resolution to the problems encountered by Beveridge 
and others in their attempts to detect the hidden periodicities in economic and 
astronomical data. Even more remarkable is the way in which the generalised 
harmonic analysis that gave rise to the concept of the spectral density of a time 
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series should prove to be wholly conformable with the alternative methods of time- 
series analysis in the time domain which arose largely as a consequence of the failure 
of the traditional methods of periodogram analysis. 

The synthesis of the two branches of time-series analysis was achieved inde- 
pendently and almost simultaneously in the early 1930s by Norbert Wiener [522] 
in America and A. Khintchine [289] in Russia. The Wiener-Khintchine theorem 
indicates that there is a one-to-one relationship between the autocovariance func- 
tion of a stationary process and its spectral density function. The relationship is 
expressed, in one direction, by writing 

^ OO 

(1.13) g{u) = — Y 7r cos(wt); 7 t = 7-r, 

T — — OQ 



where g(u) is the spectral density function and {y T ; r = 0, 1, 2, . . .} is the sequence 
of the autocovariances of the series y[t). 

The relationship is invertible in the sense that it is equally possible to express 
each of the autocovariances as a function of the spectral density: 

(1.14) y T = / cos(ujT)g(uj)duj. 

Jo 

If we set r = 0, then cos(wt) = 1, and we obtain, once more, the equation (1.12) 
which neatly expresses the way in which the variance 70 = V{y(t)} of the series 
y(t) is attributable to the constituent harmonic components; for g(cu) is simply the 
expected value of the squared amplitude of the component at frequency u>. 

We have stated the relationships of the Wiener-Khintchine theorem in terms of 
the theoretical spectral density function g(u>) and the true autocovariance function 
{7 t',t = 0,1,2,...}. An analogous relationship holds between the periodogram 
I(uJj) defined in (1.2) and the sample autocovariance function {c T ;r = 0,1,..., 
T — 1} where c T = J2(Vt ~ V){yt-T ~ y) /T- Thus, in the appendix, we demonstrate 
the identity 



T—l 

(1.15) I{iOj) = 2 c t cos(wj-t); c t = c_ T . 

t—l -T 



The upshot of the Wiener-Khintchine theorem is that many of the techniques 
of time-series analysis can, in theory, be expressed in two mathematically equivalent 
ways which may differ markedly in their conceptual qualities. 

Often, a problem which appears to be intractable from the point of view of one 
of the domains of time-series analysis becomes quite manageable when translated 
into the other domain. A good example is provided by the matter of spectral 
estimation. Given that there are difficulties in computing all T of the ordinates of 
the periodogram when the sample size is large, we are impelled to look for a method 
of spectral estimation which depends not upon smoothing the periodogram but 
upon performing some equivalent operation upon the sequence of autocovariances. 
The fact that there is a one-to-one correspondence between the spectrum and the 
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Figure 1.7. The periodogram of Wolfer’s sunspot numbers 1749-1924. 



sequence of autocovariances assures us that this equivalent operation must exist; 
though there is, of course, no guarantee that it will be easy to perform. 

In fact, the operation which we perform upon the sample autocovariances is 
simple. For, if the sequence of autocovariances {c T ; r = 0, 1, . . . , T — 1} in (1.15) is 
replaced by a modified sequence {w T c T ] r = 0, 1, . . . , T— 1} incorporating a specially 
devised set of declining weights {u> T ; r = 0, 1, . . . , T — 1}, then an effect which is 
much the same as that of smoothing the periodogram can be achieved (compare 
Figures 1.7 and 1.8). Moreover, it may be relatively straightforward to calculate 
the weighted autocovariance function. 

The task of devising appropriate sets of weights provided a major research 
topic in time-series analysis in the 1950s and early 1960s. Together with the task 
of devising equivalent procedures for smoothing the periodogram, it came to be 
known as spectral carpentry. 

The Maturing of Time-Series Analysis 

In retrospect, it seems that time-series analysis reached its maturity in the 
1970s when significant developments occurred in both of its domains. 

A major development in the frequency domain occurred when Cooley and 
Tukey [125] described an algorithm which greatly reduces the effort involved in 
computing the periodogram. The fast Fourier transform (FFT), as this algorithm 
has come to be known, allied with advances in computer technology, has enabled 
the routine analysis of extensive sets of data; and it has transformed the procedure 
of smoothing the periodogram into a practical method of spectral estimation. 

The contemporaneous developments in the time domain were influenced by 
an important book by Box and Jenkins [70]. These authors developed the time- 
domain methodology by collating some of its major themes and by applying it 
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Figure 1.8. The spectrum of the sunspot numbers calculated from 
the autocovariances using Parzen’s [383] system of weights. 



to such important functions as forecasting and control. They demonstrated how 
wide had become the scope of time-series analysis by applying it to problems as 
diverse as the forecasting of airline passenger numbers and the analysis of com- 
bustion processes in a gas furnace. They also adapted the methodology to the 
computer. 

Many of the current practitioners of time-series analysis have learnt their skills 
in recent years during a time when the subject has been expanding rapidly. Lack- 
ing a longer perspective, it is difficult for them to gauge the significance of the 
recent practical advances. One might be surprised to hear, for example, that, as 
late as 1971, Granger and Hughes [227] were capable of declaring that Beveridge’s 
calculation of the periodogram of the wheat price index (see 14.4), comprising 
300 ordinates, was the most extensive calculation of its type to date. Nowadays, 
computations of this order are performed on a routine basis using microcomputers 
containing specially designed chips which are dedicated to the purpose. 

The rapidity of the recent developments also belies the fact that time-series 
analysis has had a long history. The frequency domain of time-series analysis, to 
which the idea of the harmonic decomposition of a function is central, is an inheri- 
tance from Euler (1707-1783), d’Alembert (1717 -1783), Lagrange (1736-1813) and 
Fourier (1768-1830). The search for hidden periodicities was a dominant theme of 
nineteenth century science. It has been transmogrified through the refinements of 
Wiener’s generalised harmonic analysis [522] which has enabled us to understand 
how cyclical phenomena can arise out of the aggregation of random causes. The 
parts of time-series analysis which bear a truly twentieth-century stamp are the 
time-domain models which originate with Slutsky and Yule and the computational 
technology which renders the methods of both domains practical. 



15 



D.S.G. POLLOCK: TIME-SERIES ANALYSIS 



The effect of the revolution in digital electronic computing upon the 
practicability of time-series analysis can be gauged by inspecting the purely 
mechanical devices (such as the Henrici-Conradi and Michelson-Stratton harmonic 
analysers invented in the 1890s) which were once used, with very limited success, to 
grapple with problems which are nowadays almost routine. These devices, some of 
which are displayed in London’s Science Museum, also serve to remind us that many 
of the developments of applied mathematics which startle us with their modernity 
were foreshadowed many years ago. 

Mathematical Appendix 

Mathematical Expectations 

The mathematical expectation or the expected value of a random variable x is 
defined by 

/ OO 

xdF(x), 

-CO 

where F{x) is the probability distribution function of x. The probability distribu- 
tion function is defined by the expression F(x*) = P{x < x*} which denotes the 
probability that x assumes a value no greater than x*. If F{x) is a differentiable 
function, then we can write dF{x) = f{x)dx in equation (1.16). The function 
/(x) = dF(x)/dx is called the probability density function. 

If y(t) = {y t \ t = 0, ±1, ±2, . . .} is a stationary stochastic process, then E(y t ) = 
/i is the same value for all t. 

If S/o, - - - , Ut—i is a sample of T values generated by the process, then we may 
estimate y from the sample mean 



(1.17) 




Hvt- 



Autocovariances 

The autocovariance of lag r of the stationary stochastic process y(t) is defined 
by 

(1-18) 7 r = E{(y t - n)(y t - T - ix)}. 

The autocovariance of lag r provides a measure of the relatedness of the elements 
of the sequence y{t) which are separated by r time periods. 

The variance, which is denoted by V{y(t)} = 70 and defined by 

(1-19) 7o = E{(y t -p) 2 } 1 

is a measure of the dispersion of the elements of y{t). It is formally the autocovari- 
ance of lag zero. 
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If y t and yt- T are statistically independent, then their joint probability density 
function is the product of their individual probability density functions so that 
f(yt,yt-r) = f(yt)f(yt-r)- it follows that 

(1.20) 7r = E(y t - n)E(y t - T - y) = 0 for all r/0. 

If yo, , yr-i is a sample from the process, and if r < T, then we may estimate 
7 T from the sample autocovariance or empirical autocovariance of lag r: 

1 T_1 

( 1 - 21 ) (A = -^2(y t -y)(y t - r -y). 

t=T 



The Periodogram and the Autocovariance Function 
The periodogram is defined by 



(1.22) I(wj) = - 



T- 1 n 2 , T— 1 

Y cos (Ujt)(y t -y)[ + l Y sin ( w ^)(^ - V) \ 
t = o -I t — o -I 



The identity cos (vjt)(yt — y) = cos(u> jt)y t follows from the fact that, by 
construction, ]T) t cos(ujjt) = 0 for all j. Hence the above expression has the same 
value as the expression in (1.2). Expanding the expression in (1.22) gives 



(1.23) 



!(uj) = ^ { Y cos ( w ^) cos {vjs)(yt - y){y s - y) \ 

+ j, | sin(wjt) sin(o ; : ,s)(j/t - y)(y s - y)|, 



and, by using the identity cos (H) cos (B) + sin(vl) sin(H) = cos(H — B ), we can 
rewrite this as 

(i.24) Huj) = |;| Y^2 cos ( w A L - s Y)(yt - y)(y s -y) 

^ t s 



Next, on defining t = t — s and writing c T = ^2 t {yt — y)(yt.-T — y)/T , we can reduce 
the latter expression to 



T - 1 

(1.25) I{uj) = 2 ^ cos {cOjT )c r , 

r=l-T 



which appears in the text as equation (1.15). 
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Polynomial Methods 




CHAPTER 2 



Elements of 
Polynomial Algebra 



In mathematical terminology, a time series is properly described as a temporal 
sequence; and the term series is reserved for power series. By transforming temporal 
sequences into power series, we can make use of the methods of polynomial algebra. 
In engineering terminology, the resulting power series is described as the 2 -transform 
of the sequence. 

We shall begin this chapter by examining some of the properties of sequences 
and by defining some of the operations which may be performed upon them. Then 
we shall examine briefly the basic features of time-series models which consist of lin- 
ear relationships amongst the elements of two or more sequences. We shall quickly 
reach the opinion that, to conduct the analysis effectively, some more mathematical 
tools are needed. Amongst such tools are a variety of linear operators defined on 
the set of infinite sequences; and it transpires that the algebra of the operators 
is synonymous with the algebra of polynomials of some indeterminate argument. 
Therefore, we shall turn to the task of setting forth the requisite results from the 
algebra of polynomials. In subsequent chapters, further aspects of this algebra will 
be considered, including methods of computation. 



Sequences 

An indefinite sequence x(t ) = \x t ; t = 0, ±1, ±2, . . .} is any function mapping 
from the set of integers Z = {t = 0, ±1, ±2, . . .} onto the real line 1Z or onto the 
complex plane C. The adjectives indefinite and infinite may be used interchangeably. 
Whenever the integers represents a sequence of dates separated by a unit time 
interval, the function x(t) may be described as a time series. The value of the 
function at the point r £ Z will be denoted by x T = x(r). The functional notation 
will be used only when r G Z, which is to say when r ranges over the entire set of 
positive and negative integers. 

A finite sequence {ao, aq, . . . , a p } is one whose elements may be placed in a 
one-to-one correspondence with a finite set of consecutive integers. Such sequences 
may be specified by enumeration. Usually, the first (nonzero) element of a finite 
sequence will be given a zero subscript. A set of T observations on a time series 
x(t) will be denoted by Xo, Xi , . . . , Xt- i- Occasionally, t itself will denote a nonzero 
base index. 

It is often convenient to extend a finite sequence so that it is defined over 
the entire set of integers Z. An ordinary extension a(i ) of a finite sequence 
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{ao,ai, . . . ,a p } is obtained by extending the sequence on either side by an in- 
definite number of zeros. Thus 



(2.1) 



a(i ) 



cci, for 0 < i < p; 
0, otherwise. 



A periodic extension a{i ) of the sequence is obtained by replicating its elements 
indefinitely in successive segments. Thus 



(2.2) 



a{i) 



cti, for 0 < i < p\ 

a(imod[p+i]), otherwise, 



where (i mod [p+ 1]) is the (positive) remainder after the division of i by p+ 1. The 
ordinary extension of the finite sequence {au, aq, . . . , a p } and its periodic extension 
are connected by the following formula: 



OO 

(2.3) a(i)= Y a(i+\p+l]j). 

i=-oo 



It is helpful to name a few sequences which are especially important for analytic 
purposes. The sequence specified by the conditions 



(2.4) 



<5(t) 



1, if r = 0; 
0, if r ^ 0 



is called the unit impulse. The formulation which takes i as the index of this 
sequence and which sets 5(i — j) = 1 when i = j and S(i — j) = 0 when i ^ j 
reminds us of Kronecker’s delta. The continuous-time counterpart of the impulse 
sequence is known as Dirac’s delta. 

The unit-step sequence is specified by 



(1, if t 

( 2 - 5 ) U ( T ) = \ n ., 

10, if T 

This is related to the unit-impulse sequence by 



> 0 ; 

< 0. 

the equations 



( 2 . 6 ) 



S(t) = u(t) — u(t — 1) and 



u ( t ) = Y 6 W- 

t= — OO 



Also of fundamental importance are real and complex exponential sequences. 
A real exponential sequence is given by the function x{t) = e rt = a 4 where a = e r 
is a real number. A complex exponential sequence is given by x(t) = e lwt+ ^ with 
i = y/—l. A sinusoidal sequence is a combination of conjugate complex exponential 
sequences: 

(2.7) x(t) = pcos(ut -9)= ^{e i(u;t - e) + e"* (w ‘- e) }. 
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Here oj, which is a number of radians, is a measure of the angular velocity or angular 
frequency of the sinusoid. The parameter p represents the amplitude or maximum 
value of the sinusoid, whilst 0 is its phase displacement. 

A sequence x(t) is said to be periodic with a period of T if x(t + T) = x(t) 
for all integer values t or, equivalently, if x(t) = x(t mod T). The function x(t) = 
pcos(uit — 6) is periodic in this sense only if 2t:/u> is a rational number. If 2-k/u = T 
is an integer, then it is the period itself. In that case, its inverse / = u/2t: is the 
frequency of the function measured in cycles per unit of time. 

In some cases, it is helpful to define the energy of a sequence x(t) as the sum of 
squares of the moduli of its elements if the elements are complex valued, or simply 
as the sum of squares if the elements are real: 

(2.8) J=^|a: t | 2 . 

In many cases, the total energy will be unbounded, although we should expect it 
to be finite over a finite time interval. 

The power of a sequence is the time-average of its energy. The concept is 
meaningful only if the sequence manifests some kind of stationarity. The power of a 
constant sequence x(t) = a is just a 2 . The power of the sequence x(t) = pcos(wf) is 
ip 2 . This result can be obtained in view of the identity cos 2 {u>t) = ^{l + cos(2aA)}; 
for the average of cos(2aA) over an integral number of cycles is zero. 

In electrical engineering, the measure of the power of an alternating current is 
its so-called mean-square deviation. In statistical theory, the mean-square deviation 
of a finite sequence of values drawn from a statistical population is known as the 
sample variance. The notions of power and variance will be closely linked in this 
text. 

When the condition that 

(2.9) 5>*l<o° 

is fulfilled, the sequence x(t) = {x t ;t = 0,±1,±2, ...} is said to be absolutely 
summable. A sequence which is absolutely summable has finite energy and vice 
versa. 

There are numerous operations which may be performed upon sequences. 
Amongst the simplest of such operations are the scalar multiplication of the ele- 
ments of a sequence, the pairwise addition of the elements of two sequences bearing 
the same index and the pairwise multiplication of the same elements. Thus, if 
A is a scalar and x(t) = {x t } is a sequence, then A x(t) = {Aa: t } is the sequence 
obtained by scalar multiplication. If x{t) = {xy} and y(t) = { yt } are two se- 
quences, then x(t) + y(t) = {xt. + yt} is the sequence obtained by their addition, 
and x(t)y(t) = {x t yt} is the sequence obtained by their multiplication. 

In signal processing, a multiplication of one continuous-time signal by another 
often corresponds to a process of amplitude modulation. This entails superimposing 
the characteristics of a (continuous-time) signal y[t) onto a carrier x(t) so that in- 
formation in the signal can be transmitted by the carrier. Usually, the unmodulated 
carrier, which should contain no information of its own, has a periodic waveform. 

Also of fundamental importance are the operations linear and circular convo- 
lution which are described in the next two sections. 
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Figure 2.1. A method for finding the linear convolution of two sequences. 
The element 74 = a\0z + a 2 02 of the convolution may be formed by multiply- 
ing the adjacent elements on the two rulers and by summing their products. 



Linear Convolution 

Let {ao,ai, . . . ,a p } and {/ 3 0 , 0i, ■ ■ ■ , 0k} be two finite sequences, and consider 
forming the pairwise products of all of their elements. The products can be arrayed 
as follows: 

oto0o &o0i a o 0 2 ■ ■ ■ <+o0k 
& 1 P 0 &i0i a\0 2 ■ ■ ■ aipk 
(0 ini <+ 2 0o ot 2 0\ a 2 0 2 . . . a 2 (3k 



Otpfio tXp0 1 Op/?2 • • - Oi p 0k- 

Then a sequence 70, 71, . . . , ”/ P + q can be defined whose elements are obtained by 
summing the elements of the array along each of the diagonals which run in the 
NE-SW direction: 



70 — c*o/3o, 

71 = a oPi + a i0Ot 

(2.11) 72 = aofo + «i/5i + a 2 p 0 , 

7 p+k — Op 0 k - 

The sequence {7,} is described as the convolution of the sequences {07} and { 0 ji- 
lt will be observed that 



(2.12) 



p+k 

i 2 ' 0 = \ 


(£«/) 


( ^ 0 -j I and that 


3=0 


V j—o / 


V 0 / 


p+k 


/ P N 


. / k \ 


X]l 7 ll< | 
3=0 


1 Ek-i 

v j= 0 > 


)(£l&l)- 

v i= 0 7 
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Example 2.1. The process of linear convolution can be illustrated with a simple 
physical model. Imagine two rulers with adjacent edges (see Figure 2 . 1 ). The lower 
edge of one ruler is marked with the elements of the sequence {ao, ai, . . . , a p } at 
equally spaced intervals. The upper edge of the other ruler is marked with the 
elements of the reversed sequence {Pk, ■ ■ ■ , Pi, Po} with the same spacing. At first, 
the rulers are placed so that ag is above Pq. The pair (ao,/ 3 o) is written down 
and the product 70 = «o/?o is formed. Next the lower ruler is shifted one space to 
the right and the pairs (ao,/3i) and (ai,/?o) are recorded from which the sum of 
products 71 = aoPi + ai/3o is formed. The lower ruler is shifted to the right again 
and 72 is formed. The process continues until the final product J p +k = ot p Pk is 
formed from the pair ( a p ,Pk ). 

The need to form the linear convolution of two finite sequences arises very 
frequently, and a simple procedure is required which will perform the task. The 
generic element of the convolution ol {ao,a\, ... ,a p } and {Po, Pi, ■ • • , Pk}, is given 
by 



S 

= T aiPj-i, where 
( 2 . 13 ) U 

r = max(0,j — k) and s = min(p, j). 



Here the restriction r < i < s upon the index of the summation arises from the 
restrictions that 0 < i < p and that 0 < (j — i) < k which apply to the indices of 
ai and Pj-i- 

In the following procedure, which implements this formula, the elements 7 j are 
generated in the order of j = p + k, . . . , 0 and they are written into the array beta 
which has hitherto contained the elements of {Pj}. 



(2.14) procedure Convolution^ ar alpha, beta : vector ; 

p , k : integer ); 



var 

gamma : real ; 
i,j, r, s : integer ; 

begin 

for j : = p + k downto 0 do 
begin {j} 

s := Min(j,p); 
r := Max(0,j — k)-, 
gamma := 0.0; 

for * := r to s do 

gamma := gamma + alpha[i\ * beta[j — i\; 
beta[j] := gamma ; 
end: {;/'} 

end; {Convolution} 
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Some care must be exercised in extending the operation of linear convolution 
to the case of indefinite sequences, since certain conditions have to be imposed to 
ensure that the elements of the product sequence will be bounded. The simplest 
case concerns the convolution of two sequences which are absolutely summable: 

(2.15) If a(i) = {a,} and /3(i) = {pi} are absolutely summable sequences 

such that < oo and Y} |/3.;| < oo, then their convolution product, 

which is defined by 



i— — oo i=— oo 

is also an absolutely summable sequence. 

Here the absolute summability of the product sequence, which entails its bounded- 
ness, can be demonstrated by adapting the inequality under (2.12) to the case of 
infinite sequences. 

Circular Convolution 

Indefinite sequences which are obtained from the periodic extension of finite 
sequences cannot fulfil the condition of absolute summability; and the operation of 
linear convolution is undefined for them. However, it may be useful, in such cases, 
to define an alternative operation of circular convolution. 

(2.16) Let a(i) = {ay} and /3(i) = {Pi} be the indefinite sequences 
which are formed by the periodic extension of the finite sequences 
{ao> Q!i, ... , a n -i} and {/3 q, P i, • • • , Pn- 1} respectively. Then the cir- 
cular convolution of a(i) and f3{i) is a periodic sequence defined by 

n— 1 n— 1 

7 (j) = “^0 -i) = ^2 PMi - *)■ 

2—0 2—0 

To reveal the implications of the definition, consider the linear convolution of the 
finite sequences {ay} and {Pj} which is a sequence {70, 71, 7271-2}- Also, let 

oij = oi(jmoAn) and (3j = /3(y m odn) denote elements of a{i) and /3(i). Then the 
generic element of the sequence 7 (J) is 

j n — 1 

7 j='52°‘iPj-i+ a A-i 

2—0 2=J + 1 

(2.17) J "- 1 

^ ' C^iPj — i T ^ ' ®iPj-\-n—i 

i— 0 i—j+1 

= To+Tj+n- 

The second equality depends upon the conditions that ciy = a^ when 0 < * < n, 
that Pj-i = Pj-i when 0 < (j — i) < n and that Pj-i = P(j~i) mo An — Pj+n-i when 
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Figure 2.2. A device for finding the circular convolution of two sequences. 
The upper disc is rotated clockwise through successive angles of 30 degrees. 
Adjacent numbers on the two discs are multiplied and the products are 
summed to obtain the coefficients of the convolution. 



— n < (j — i) < 0. Thus it can be seen that j(j) represents the periodic extension 
of the finite sequence {70, 71 , . . . , 77,-1} wherein 

(2.18) 7j = 7j + 7j+n for j = 0,...,n-2 and 7„_i = 7„_i. 



Example 2.2. There is a simple analogy for the process of circular convolution 
which also serves to explain the terminology. One can imagine two discs placed one 
above the other on a common axis with the rim of the lower disc protruding (see Fig- 
ure 2.2). On this rim, are written the elements of the sequence {ao, a 1, ■ ■ ■ , a n -i} at 
equally spaced intervals in clockwise order. On the rim of the upper disc are written 
the elements of the sequence {/3q,/3i, . . . ,(3 n - 1} equally spaced in an anticlockwise 
order. The circular disposition of the sequences corresponds to the periodic nature 
of the functions a(i) and (3{i ) defined in (2.16). 

At the start of the process of circular convolution, ao and /?o are in alignment, 
and the pairs (ao,/3o)> (aq, /3 n _i), ■ • ■ , (a n -i, /?i) are read from the discs. Then, 
the upper disc is turned clockwise through an angle 2n/n radians and the pairs 
(ao, /3i), (oq, /3q), . . . , (a n _i, P 2 ) are read and recorded. The process continues until 
the (n — l)th turn when the pairs (ao, (<ai,/3 n _2), • • ■ , (a n -i,Po) are read. 

One more turn would bring the disc back to the starting position. From what has 
been recorded, one can form the products 70 = a o/?o + ®if3n-i + • • • + a n _i/3i, 
71 = ao/3i + aiflo + • • • + a n -\f32, ■ ■ ■ ,7n-i = <ao/3n-i + o-iPn-2 + • • • + a ra _i/?o 
which are the coefficients of the convolution. 

The Pascal procedure which effects the circular convolution of two sequences 
is a straightforward one: 

(2.19) procedure Circonvolve(alpha,betci : vector ; 

var gamma : vector ; 
n : integer ); 
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var 

i,j , k : integer ; 



begin 

for j := 0 to n — 1 do 
begin {j} 

gamma[j ] := 0.0; 

for i := 0 to n — 1 do 
begin 
k := j - i\ 
if k < 0 then 
k := k + n; 

gamma[j ] := gamma[j ] + alpha[i\ * beta[k]\ 

end; 
end; {j} 

end; {deconvolve} 



Time-Series Models 

A time-series model is one which postulates a relationship amongst a number 
of temporal sequences or time series. Consider, for example, the regression model 

(2.20) y(t) = Px(t)+e(t), 

where x(t) and y(t ) are observable sequences indexed by the time subscript t and 
e(f) is an unobservable sequence of independently and identically distributed ran- 
dom variables which are also uncorrelated with the elements of the explanatory 
sequence of x(t). The purely random sequence e{t) is often described as white 
noise. 

A more general model is one which postulates a relationship comprising any 
number of consecutive elements of x(t), y(t) and e{t). Such a relationship is ex- 
pressed by the equation 

p k q 

(2.21) ^2 aiy{t - i) = ^2 & iX (* _ *) + ^2 Ti£(t - i), 

2—0 2 — 0 2 = 0 

wherein the restriction a 0 = 1 is imposed in order to identify y{t ) as the output of 
the model. The effect of the remaining terms on the LHS is described as feedback. 
Any of the sums in this equation can be infinite; but, if the model is to be viable, 
the sequences of coefficients {«i}, {/3i} and {yq} must depend on a strictly limited 
number of underlying parameters. Notice that each of the terms of the equation 
represents a convolution product. 

A model which includes an observable explanatory sequence or signal sequence 
x(t) is described as a regression model. When x(t ) is deleted, the simpler uncondi- 
tional linear stochastic models are obtained. Thus the equation 

p q 

(2.22) ^2 a iy(t ~ i) = ^ lR£{t - i) 

2=0 2=0 
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represents a so-called autoregressive moving-average (AR.MA) process. When a* = 
0 for all i > 0, this becomes a pure moving-average (MA) process. When /i, = 0 
for all i > 0, it becomes a pure autoregressive (AR) process. 

Transfer Functions 

Temporal regression models are more easily intelligible if they can be repre- 
sented by equations in the form of 

(2.23) y(t) = ^2 - z) + i£(t - i), 

£>0 i >0 

where there is no lag scheme affecting the output sequence y(t). This equation 
depicts y(t) as a sum of a systematic component hit) = Y u>ix(t—i) and a stochastic 
component 77 (f) = Y^P^i t ~~ *)• Both of these components comprise transfer- 
function relationships whereby the input sequences x{t) and e(f) are translated, 
respectively, into output sequences h(t) and 77 (f). 

In the case of the systematic component, the transfer function describes how 
the signal x(f) is commuted into the sequence of systematic values which explain a 
major part of y{t) and which may be used in forecasting it. 

In the case of the stochastic component, the transfer function describes how a 
white-noise process e(f), comprising a sequence of independent random elements, is 
transformed into a sequence of serially correlated disturbances. In fact, the elements 
of h{t) represent efficient predictors of the corresponding elements of y{t) only when 
77 (f) = t/>o e(f) is white noise. 

A fruitful way of characterising a transfer function is to determine the response, 
in terms of its output, to a variety of standardised input signals. Examples of such 
signals, which have already been presented, are the unit-impulse <5 (f) , the unit-step 
u(t) and the sinusoidal and complex exponential sequences defined over a range of 
frequencies. 

The impulse response of the systematic transfer function is given by the se- 
quence h(t) = Yi UiS(t — «)• Since i € {0, 1,2,.. .}, it follows that h(t) = 0 for all 
f < 0. By setting t = {0, 1,2,.. .}, a sequence is generated beginning with 

ho = u> 0 , 

(2.24) hi = cui , 

I12 = U >2 ■ 

The impulse-response function is nothing but the sequence of coefficients which 
define the transfer function. 

The response of the transfer function to the unit-step sequence is given by 
h(t) = Yi u i u {t ~ *)• By setting t = { 0 , 1,2 , . . .}, a sequence is generated which 
begins with 

ho = ijJq ) 

(2.25) hi = u>q + u>i , 

h. 2 = U> 0 + iOl + U)2- 

Thus the step response is obtained simply by cumulating the impulse response. 
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In most applications, the output sequence h(t) of the transfer function should 
be bounded in absolute value whenever the input sequence x(t) is bounded. 
This is described as the condition of bounded input- bounded output (BIBO) 
stability. 

If the coefficients {wq; uj\, • ■ ■ > w p } of the transfer function form a finite se- 
quence, then a necessary and sufficient condition for BIBO stability is that \uii < oo 
for all i, which is to say that the impulse-response function must be bounded. 
If . . .} is an indefinite sequence, then it is necessary, in addition, that 

< oo, which is the condition that the step-response function is bounded. 
Together, the two conditions are equivalent to the single condition that Y \uJi\ < oo, 
which is to say that the impulse response is absolutely summable. 

To confirm that the latter is a sufficient condition for stability, let us consider 
any input sequence x(t) which is bounded such that |a;(t)| < M for some finite M. 
Then 



(2.26) 



\h(t)\ = Q2(AJix(t - i) 




< oo, 



and so the output sequence h(t) is bounded. To show that the condition is necessary, 
imagine that Y l w i| is unbounded. Then a bounded input sequence can be found 
which gives rise to an unbounded output sequence. One such input sequence is 
specified by 



(2.27) 



This gives 
(2.28) 



X-i 



UJi 

I u>i 

0 , 



if ^ ^ 0; 

if UJi = 0. 



h 0 = ^2 ujiX-i 



i^i> 



and so h(t) is unbounded. 

A summary of this result may be given which makes no reference to the specific 
context in which it has arisen: 



(2.29) The convolution product h{t) = Y u i x {t — *), which comprises a 
bounded sequence x(t) = {xt}, is itself bounded if and only if the 
sequence {uJi} is absolutely summable such that Yi l w i| < °°- 

In order to investigate the transfer-function characteristics of a relationship in 
the form of the general temporal model of equation (2.21), it is best to eliminate 
the lagged values of the output sequence y(t) which represent feedback. This may 
be done in a number of ways, including a process of repeated substitution. 

A simple example is provided by the equation 

(2.30) y(t) = (jjy{t - 1) + j3x(t) + e(t). 



32 



2: ELEMENTS OF POLYNOMIAL ALGEBRA 



A process of repeated substitution gives 

y{t) = - 1) + 0x(t ) + e(t) 

= <f) 2 y(t - 2) + /3{x(t) + - 1)} + e(t) + - 1) 

(2.31) : 

= 4> n y{t — n) + 0{x(t) + 4>x(t — 1) + • • • + (j) n 1 x(t — n + 1)} 
+ s(t) + (f>s{t — !) + ••■ + (j> n 1 s(t — n + 1). 



If \(f>\ < 1, then lim(n — » oo)<p n = 0; and it follows that, if x(t) and e(t) are bounded 
sequences, then, as the number of repeated substitutions increases indefinitely, the 
equation will tend to the limiting form of 

OO OO 

(2.32) y(t) = 0^2 <l> l x(t - i) + ^ 4>' e(t - i), 

i = 0 i—O 

which is an instance of the equation under (2.23). 

For models more complicated than the present one, the method of repeated 
substitution, if pursued directly, becomes intractable. Thus we are motivated to 
use more powerful algebraic methods to effect the transformation of the equation. 



The Lag Operator 

The pursuit of more powerful methods of analysis begins with the recognition 
that the set of all time series {x(t): : t G Z,x £ 1Z} represents a vector space. Various 
linear transformations or operators may be defined over the space. The lag operator 
L , which is the primitive operator, is defined by 

(2.33) Lx(t) = x(t- 1). 

Now, L{Lx(t)} = Lx(t — 1) = x(t — 2); so it makes sense to define L 2 by L 2 x(t) = 
x(t — 2). More generally, L k x{t) = x(t — k ) and, likewise, L~ k x(t ) = x(t + k). 
Other important operators are the identity operator / = L , the annihilator or zero 
operator 0 = 1 — 1, the forward-difference operator A = L^ 1 — I, the backwards- 
difference operator V = LA = I — L and the summation operator S = (I + L + 
L 2 + •••)• 

The backwards-difference operator has the effect that 

(2.34) Vx(f) = x(t) — x(t — 1), 
whilst the summation operator has the effect that 

OO 

(2.35) Sx(t) = x(t — i ). 

i=0 

These two operators bear an inverse relationship to each other. On the one 
hand, there is the following subtraction: 

s = i + l + l 2 + 

(2.36) £5* = L + L 2 + • • • 

S — LS = I. 
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This gives S(I — L) = , SV = /, from which S = V -1 . The result is familiar from 
the way in which the sum is obtained of a convergent geometric progression. On 
the other hand is expansion of S = 1/(1 — L). A process of repeated substitution 
gives rise to 



S=I + LS 

(2.37) =I + L + L 2 S 

= 1 + L + L 2 + L 3 S. 

If this process is continued indefinitely, then the original definition of the summation 
operator is derived. The process of repeated substitution is already familiar from 
the way in which the equation under (2.30), which stands for a simple temporal 
regression model, has been converted to its transfer-function form under (2.32). 

Another way of expanding the operator S = 1/(1 — L) is to use the algorithm 
of long division: 

(2.38) I + L + L 2 + --- 

i -lJI 

I -L 
L 

L-L 2 



L 2 

L 2 -L 3 

If this process is stopped at any stage, then the results are the same as those from 
the corresponding stage of the process under (2.37). The binomial theorem can also 
be used in expanding S' = (/ — L) _1 . 

To all appearances, the algebra of the lag operator is synonymous with ordinary 
polynomial algebra. In general, a polynomial of the lag operator of the form p(L) = 
Po + piL + ■ ■ ■ + p n L n = PiL 1 has the effect that 

p(L)x(t) =pox(t) + pix(t — 1) + • • • + p n x(t — n) 

(2.39) -A 

= - *)• 
i=0 

The polynomial operator can be used to re-express the temporal regression 
model of (2.21) as 



(2.40) a(L)y(t) = (3(L)x(t) + p(L)e(t). 

In these terms, the conversion of the model to the transfer-function form of equation 
(2.23) is a matter of expanding the rational polynomial operators f3(L)/a(L) and 
p(L)/a(L) in the expression 



(2.41) 



,(t) = ffiko + £(t) . 



a(L)' 



a(L) 
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We shall be assisted in such matters by having an account of the relevant algebra 
and of the corresponding algorithms readily to hand. 

Algebraic Polynomials 

A polynomial of the ptlr degree in a variable z, which may stand for a real or a 
complex-valued variable, or which may be regarded as an indeterminate algebraic 
symbol, is an expression in the form of 

(2.42) 01 ( 2 ) = a o T ot\Z + • • • + cx p z p , 

where it is understood that ao,a p ^ 0. When z £ C is a complex- valued variable, 
a(z) may be described as the ^-transform of the sequence {ao,ai, . . . ,ai}. From 
another point of view, a(z) is regarded as the generating function of the sequence. 

Let a(z) = ao + a\z + ■ ■ ■ + a p z p and (3{z) = p 0 + f3\z + • • • + PkZ k be two 
polynomials of degrees p and k respectively. Then, if k > p, their sum is defined by 

a(z) + (3(z) = (a 0 + po) + (aq + f3i)z + • • • + ( a p + P p )z p 

[ ’ + Pp+lZ P+1 + h PkZ k ■ 

A similar definition applies when k < p. 

The product of the polynomials a{z) and (3{z) is defined by 

a( z )P(z) = aoPo + (&oPi + &iPo) z + • • • + ct p (3k zP+k 

(2.44) = 7o + H z + 72 2 2 H f 7 P +k zP+k 

= 7 {z)- 

The sequence of coefficients { 7 .;} in the product is just the convolution of the se- 
quences {cti} and {Pi} of the coefficients belonging to its factors. 

These operations of polynomial addition and multiplication obey the simple 
rules of (i) associativity, (ii) commutativity and (iii) distributivity which are found 
in arithmetic. If a(z), P{z) and 7 ( 2 ) are any three polynomials, then 

(i) {a{z)P{z)}'y{z) = a{z){P{z)'y{z)}, 

{a(z) + P(z)} + 7 (z) = a(z) + {P(z) + 7 (z)}; 

(2.45) (ii) a(z) + P(z) = P(z) + a(z), 

a(z)P(z) = P(z)a(z); 

(iii) a(z){P(z) + 7 ( 2 )} = a(z)P(z) + a(z)j(z). 

Periodic Polynomials and Circular Convolution 

If the polynomial argument z° is a periodic function of the index j, then the set 
of polynomials is closed under the operation of polynomial multiplication. To be 

precise, let a(z) = a 0 +aiZ-\ ba n _i.z ” -1 and P(z) = P 0 + piz~\ b/^n-i-Z ™ -1 

be polynomials of degree n — 1 at most in an argument which is n-periodic such 
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that z 3+n = z 3 for all j, or equivalently z T j = z T (j mod n ). Then the product 
of 7 {z) = a(z)/ 3 (z) is given by 

(2 46 ^ l( z ) = 7 o + 7 iz + ■ ■ ■ + l 2 n- 2 z 2n ~ 2 

= 7 o + 7 i z-\ hVi^ 1 ! 

where the elements of {70, 71, ... , 7271-2} are the products of the linear convolution 
of {a 0 , at \, . . . , a n -i} and {Po,Pi, • • • ,Pn- 1}, and where 

( 2 . 47 ) 7 j = 7 j + 7 i+n for j = 0 , ...,n -2 and 7„_i = 7„_i. 

These coefficients (70, 7 i, ■ ■ • , 77.- 1} may be generated by applying a pro- 
cess of circular convolution to sequences which are the periodic extensions of 
{aiO; aUj • ■ • 1 OL n - 1} and {/?o> Pi, ■ ■ ■ , Pn-i}- 

The circular convolution of (the periodic extensions of) two sequences 
{ao, «i, . . . , a n _ 1} and {Pq, Pi, ■ ■ ■ , Pn- 1} may be effected indirectly via a method 
which involves finding their discrete Fourier transforms. The Fourier transforms 
of the sequences are obtained by evaluating their ^-transforms a(z), P(z) at n 
distinct points which are equally spaced around the circumference of a unit circle 
in the complex plane. These points {zk',k = 0 , ...,n— 1 } come from setting 
Zk = exp(— i 2 irk/n). From the values { a(zk)', k = 0 , 1 } and { p{zk)', 

k = 0, . . . , n — 1}, the corresponding values {7 (zk) = a(zk)P(zk)', k = 0, . . . , n — 1} 
can be found by simple multiplication. The latter represent n ordinates of the poly- 
nomial product whose n coefficients are being sought. Therefore, the coefficients 
can be recovered by an application of an (inverse) Fourier transform. 

It will be demonstrated in a later chapter that there are some highly effi- 
cient algorithms for computing the discrete Fourier transform of a finite sequence. 
Therefore, it is practical to consider effecting the circular convolution of two se- 
quences first by computing their discrete Fourier transforms, then by multiplying 
the transforms and finally by applying the inverse Fourier transform to recover the 
coefficients of the convolution. 

The Fourier method may also be used to affect the linear convolution of two 
sequences. Consider the sequences { a 0, aq, . . . , a p } and {p 0 , Pi,..., PP\ whose ,2- 
transforms may be denoted by a(z) and P(z). If z 3 is periodic in j with a period 
of n > p + k, then 

( 2 . 48 ) a(z)P(z) = 70 + 71 z -\ f 7 p+k z p+k 

resembles the product of two polynomials of a non-periodic argument, and its co- 
efficients are exactly those which would be generated by a linear convolution of the 
sequences. The reason is that the degree p+k of the product is less than the period 
n of the argument. 

In the context of the discrete Fourier transform, the period of the argument 2 
corresponds to the length n of the sequence which is subject to the transformation. 
In order to increase the period, the usual expedient is to extend the length of the 
sequence by appending a number of zeros to the end of it. This is described as 
“padding” the sequence. 
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Consider the padded sequences {ckoi a i, ■ ■ ■ , o ra _i} and {/3 0 , Pi, ■ ■ ■ , P n - 1} i n 
which a p +\ = ■ ■ ■ = o n _i = 0 and Pk+i = • ■ ■ = p n -i = 0. Let their z-transforms 
be denoted by d(z) and p{z), and let the period of z be n > p+k. Then the product 
7 (z) = a(z)/3(z), will entail the coefficients 70 = 70, 7i = 7i> • •• ,7 p+k = 7 p+k of 
the linear convolution of {«o, a \, . . . , a p } and {Po, Pi, ■ ■ ■ , Pk} together with some 
higher-order coefficients which are zeros. Nevertheless, these are the coefficients 
which would result from applying the process of circular convolution to the padded 
sequences; and, moreover, they can be obtained via the Fourier method. 

The result can be made intuitively clear by thinking in terms of the physical 
model of circular convolution illustrated in Figure 2.2. If the sequences which are 
written on the rims of the two discs are padded with a sufficient number of zeros, 
then one sequence cannot engage both the head and the tail of the other sequence 
at the same time, and the result is a linear convolution. 



Polynomial Factorisation 

Consider the equation a.^ + a\Z + a^z 1 = 0. This can be factorised as a 2 (z — 
Ai)(z — A2) where Ai, A2 are the roots of the equation which are given by the 
formula 



(2.49) 



A = 



-ai ± \Jol\ - 4a 2 ao 
2a 2 



If a\ > 4a 2 a:o, then the roots Ai, A 2 are real. If af = 4a 2 ao, then Ai = A 2 . If af < 
4a 2 ao, then the roots are the conjugate complex numbers A = a + ip, A* = a — ip 
where i = 

It is helpful to have at hand a Pascal procedure for finding the roots of a 
quadratic equation: 



(2.50) procedure QuadraticRoots(a, b, c : real); 

var 

discriminant, rootl,root2, modulus : real; 



begin 

discriminant := Sqr(b) — 4 * a * c; 
if ( discriminant > 0) and (a <> 0) then 
begin 

roofl := (— b+ Sqrt(discriminant )) / (2 * a); 
root2 := (— b — Sqrt(discriminant )) / (2 * a); 

Writeln(' Root( 1) = ' ,root\ : 10 : 5); 

Writeln(' Root{ 2) = ' ,root2 : 10 : 5); 

end; 

if ( discriminant = 0) and (a <> 0) then 
begin 

roofl := —b / (2 * a); 

Writeln('The roots coincide at the value = ' ,rootl : 10 : 5); 

end; 

if ( discriminant < 0) and (a <> 0) then 
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begin 

root 1 := —6/(2 * a); 
root2 := Sqrt(— discriminant) / (2 * a); 
modulus := Sqrt(Sqr(rootl) + Sqr(root2))\ 
Writeln('We have conjugate complex roots ' ); 
Writeln('The real part is ' ,rootl : 10 : 5); 
Writeln('The imaginary part is ' ,root2 : 10 : 5); 
Writeln('The modulus is ', modulus : 10 : 5); 
end; 

end; {Quadratic Roots} 



Complex Roots 

There are three ways of representing the conjugate complex numbers A and A*: 



The three representations are called, respectively, the Cartesian form, the trigono- 
metrical form and the exponential form. The parameter p = |A| is the modulus of 
the roots and the parameter 9, which is sometimes denoted by 9 = arg(A), is the 
argument of the exponential form. This is the angle, measured in radians, which A 
makes with the positive real axis when it is interpreted as a vector bound to the 
origin. Observe that 9 is not uniquely determined by this definition, since the value 
of the tangent is unaffected if 2mr is added to or subtracted from 9 , where n is 
any integer. The principal value of arg(A), denoted Arg(A), is the unique value of 
9 £ (—7 r, 7r] which satisfies the definition. 

The Cartesian and trigonometrical representations are understood by consid- 
ering the Argand diagram (see Figure 2.3). The exponential form is understood by 
considering the series expansions of cosf? and ism 9 about the point 9 = 0 : 



(2.51) 



A = a + i/3 = p { cos 9 + i sin 9) = pe 16 , 
A* = a — i/3 = p{ cos 9 — i sin 9) = pe~ l6 



Here there are 



(2.52) 



p = \J a 2 + /3 2 and tan 9 = (3/a. 



(2.53) 





Adding the series gives 



(2.54) 




Likewise, subtraction gives 
(2.55) 



cos 9 — ism 9 = e 
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Figure 2.3. The Argand diagram showing a complex 
number X = a + i/3 and its conjugate A* = a — ip. 



The equations (2.54) and (2.55) are known as Euler’s formulae. The inverse formu- 
lae are 

(2.56) 
and 

(2.57) sin 0 = ~-(e ie - e ~ ie ) = . 

In some computer languages — for example, in FORTRAN — complex numbers 
correspond to a predefined type; and the usual operations of complex arithmetic 
are also provided. This is not the case in Pascal; and it is helpful to define the 
complex type and to create a small library of complex operations. It is convenient 
to use a record type: 

(2.58) type 

complex = record 
re, im : real ; 

end; 

The modulus of a complex number defined in (2.52) and its square are provided 
by the following two functions: 

(2.59) function Cmod{a : complex ) : real ; 

begin 

Cmod := Sqrt(Sqr(a.re) + Sqr(a.im))' 
end; {Cmod} 
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function Cmodsqr(a : complex ) : real, 

begin 

Cmodsqr := Sqr(a.re) + Sqr(a.im); 
end; {Cmodsqr} 

The addition of a pair of complex numbers is a matter of adding their real and 
imaginary components: 

(2.60) (a + i/3) + (7 + iS) = (a + 7) + i{(3 + 5). 

Functions are provided both for addition and for subtraction: 

(2.61) function Cadd(a, b : complex ) : complex', 

var 

c : complex', 

begin 

c.re := a.re + b.re; 
c.irn := a.irn + b.im; 

Cadd := c; 
end; {Cadd} 

function Csubtract(a, b : complex) : complex; 

var 

c : complex; 

begin 

c.re := a.re — b.re; 
c.im := a.im — b.im. ; 

C subtract := c; 
end; {Csubtract} 

The product of the numbers 

A = a + i/3 = p ( cos 9 + isin60 = pe , 

(2.62) . 
p = 7 + id = k(cos w + i sin uj) = ne lu} 



is given by 



(2.63) 



Xp = 07 — pd + i{aS + pj) 

= pn{( cos 9 cos u> — sin 9 sin u>) + i ( cos 9 sin to + sin 9 cos to ) } 
= pn{ cos {9 + u>) + i sin(0 + u >) } 



= P Ke i(e+uj) 



where two trigonometrical identities have been used to obtain the third equality. 

In the exponential form, the product of the complex numbers comprises the 
product of their moduli and the sum of their arguments. The exponential repre- 
sentation clarifies a fundamental identity known as DeMoivre’s theorem: 

(2.64) {p(cos0 + isin#)}" = p n { cos {n0) + «sin(n0)}. 
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In exponential form, this becomes {pe l8 } n = p n e m8 . 

For the purposes of computing a complex multiplication, the Cartesian repre- 
sentation is adopted: 

(2.65) function Cmultiply(a , b : complex) : complex ; 

var 

c : complex', 

begin 

c.re := a.re * b. re — a.im * b.im ; 
c.im := a.im * b. re + b.im. * a.re; 

Cmultiply := c; 
end; {Cmultiply} 



The inverse of the number a + i (3 is 



( 2 . 66 ) 



(a + i{3) 1 



a — i(3 
a 2 + f3 2 ' 



This is obtained from the identity A -1 = A*/(A*A). A formula for the division of 
one complex number by another follows immediately; and the trigonometrical and 
polar forms of these identities are easily obtained. Separate code is provided for 
the operations of inversion and division: 



(2.67) function Cinverse(a : complex) : complex ; 

var 

c : complex', 

begin 

c.re := a.re/(Sqr(a.re) + Sqr(a.im)); 
c.im := — a.im /{Sqr {a.re ) + Sqr(a.im)); 

Cinverse := c; 
end; { Cinverse } 

function Cdivide(a, b : complex) : complex ; 

var 

c : complex', 

begin 

c.re := (a.re * b. re + a.im. * b.im) / {Sqrlb.re) + Sqr(b.im)); 
c.im := ( a.im * b. re — b.im * a.re) / {Sqrlb.re) + Sqr(b.im))', 

Cdivide := c; 
end; {Cdivide} 

Finally, there is the code for computing the square root of a complex number. 
In this case, the polar representation is used: 

(2.68) function Csqrt(a : complex) : complex', 

const 
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virtualZero = IE — 12; 
pi = 3.1415926; 



var 

rho, theta : real ; 
c : complex', 



begin {complex square root} 

rho := Sqrt(Sqr(a.re) + Sqr(a.im))', 
if Abs(a.re) < virtualZero then 
begin 

if a.im < 0 then 

theta := pi / 2 

else 

theta := —pi/2 

end 

else if a.re < 0 then 

theta := ArcTan[a.im/ a.re) + pi 
else 

theta := Arctan(a.im / a.re)] 
c.re := Sqrt(rho) * Cos{theta/ 2); 
c.im := Sqrt(rho ) * Sinltheta/ 2); 
Csqrt := c; 

end; {Csgrf : complex square root} 



The Roots of Unity 

Consider the equation z n = 1. This is always satisfied by z = 1; and, if n is 
even, then z = — 1 is also a solution. There are no other solutions amongst the set 
of real numbers. The solution set is enlarged if z is allowed to be a complex number. 
Let z = pe lB = p{cos(9) + isin(0)}. Then z n = p n e ien = 1 implies that p n = 1 and 
therefore p = 1, since the equality of two complex numbers implies the equality 
of their moduli. Now consider z n = e %Bn = cos {nO) + i sin(n$). Equating the real 
parts of the equation z n = 1 shows that cos (nd) = 1 which implies that nd = 2irk, 
where k is any integer. Equating the imaginary parts shows that sin(n0) = 0 which, 
again, implies that n8 = 2nk. Therefore, the solutions take the form of 



(2.69) 




2irk 2irk 

= cos b i sm . 



n n 



Such solutions are called roots of unity. 

Since cos{27rfc/n} and sin{27rA;/n} are periodic functions with a period of k = 
n, it makes sense to consider only solutions with values of k less than n. Also, if z 
is a root of unity, then so too is z* . The nth roots of unity may be represented by 
n equally spaced points around the circumference of the unit circle in the complex 
plane (see Figure 2.4). 

The roots of unity are entailed in the process of finding the discrete Fourier 
transform of a finite sequence. Later in the present chapter, and in Chapter 14, 
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Figure 2.4. The 6th roots of unity inscribed in the unit circle. 



where we consider the discrete Fourier transform of a sequence of data points y t ;t = 
0, . . . , T — 1, we adopt the notation 

(2.70) W ]! = exp : t = 0,...,T-l 

to describe the T points on the unit circle at which the argument z° is evaluated. 
The Polynomial of Degree n 

Now consider the general equation of the ?rtlr degree: 

(2.71) ol o + ot\z + • ■ • + a n z n = 0. 

On dividing by a n , a monic polynomial is obtained which has a unit associated 
with the highest power of z. This can be factorised as 

(2.72) (j» — Ai)(« — A 2 ) - - - (« — A„) = 0, 

where some of the roots may be real and others may be complex. If the coefficients 
of the polynomial are real-valued, then the complex roots must come in conjugate 
pairs. Thus, if A = a + i/3 is a complex root, then there is a corresponding root 
A * = a — i/3 such that the product (z — A) (2 — A*) = z 2 + 2az+ (a 2 + /3 2 ) is real and 
quadratic. When the n factors are multiplied together, we obtain the expansion 

(2.73) 0 = 2 "-^ XiZ n ~ 1 + Y, KXjZ n ~ 2 - • • • (-l) n (A 1 A 2 • • • A n ). 

i i^j 

This can be compared with the expression (ao/a n ) + ( a\/a n )z + • • ■ + z n = 0. By 
equating coefficients of the two expressions, it is found that ( a 0 /a n ) = (— l) n A , 
or, equivalently, 



(2.74) a n = a 0 [](-A I )- 1 - 

i—1 
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Thus the polynomial may be expressed in any of the following forms: 



(2.75) 



^2,aiZ l = a n ~^\{z - \i) 

=a 0 n(-^r i n(*-*o 

= ao na -*/*). 



The following procedure provides the means for compounding the coefficients 
of a monic polynomial from its roots. The n roots are contained in a complex array 
lambda. When all the complex roots are in conjugate pairs, the coefficients become 
the real elements of the complex array alpha. 

(2.76) procedure RootsToCoe f ficients(n : integer, 

var alpha , lambda : complexVector ) ; 



var 

j, k : integer, 
store : complex', 

begin {RootsToCoe f ficients} 
alpha[0].re := 1.0; 
alpha[0\.im := 0.0; 

for k := 1 to n do 
begin {k} 

alpha[k\.im := 0.0; 
alpha[k].re := 0.0; 

for j := k downto 1 do 
begin {j} 

store := Cmultiply(lambda[k\, alpha[j]); 
alpha[j } := C subtract (alpha[j — l],sto?'e); 
end; {j} 

alpha[ 0] := Cmultiply(lambda[k\,alpha[ 0]); 
alpha[0\.re := — alpha[0\.re', 
alpha[0\.im := — alpha[0].im 
end; {k} 

end; {RootsToCoe f ficients} 

Occasionally it is useful to consider, instead of the polynomial a(z ) of (2.71), 
a polynomial in the form 

(2.77) Q ;/ (' 2 ) = u oz n + ot\z n 1 + • • • + OL n —\z + a n . 

This has the same coefficients as a(z), but is has declining powers of z instead of 
rising powers. Reference to (2.71) shows that a'(z) = z n a(z~ 1 ). 
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If A is a root of the equation a(z) = X) a i zl = 0, then p = 1/A is a root of the 
equation a'(z) = = 0. This follows since Y^ a ih n ~ l = A 4 ” S a ih~ l = 0 

implies that cnipT 1 = )T) o^A* = 0. Confusion can arise from not knowing which 
of the two equations one is dealing with. 

Another possibility, which may give rise to confusion, is to write the factorisa- 
tion of a{z) in terms of the inverse values of its roots which are the roots of a(z _1 ). 
Thus, in place of the final expression under (2.75), one may write 

(2.78) y ^aiZ 1 = aoTT(l ~ A^~)- 

Since it is often convenient to specify a(z ) in this manner, a procedure is provided 
for compounding the coefficients from the inverse roots. In the procedure, it is 
assumed that ao = 1. 

(2.79) procedure Inver seRootsToCoef fs(n : integer ; 

var alpha, mu : complexVector ); 



var 

j, k : integer; 
store : complex; 

begin 

alpha[0\.re := 1.0; 
alpha[0\.im := 0.0; 

for k := 1 to n do 
begin {k} 

alpha[k\.im := 0.0; 
alpha[k\.re := 0.0; 

for j := k downto 1 do 
begin {j} 

store := Cmultiply(mu[k\, alpha[j — 1]); 
alpha[j } := Csubtract(alpha[j], store); 
end ;{j} 
end; {fc} 

end ; { Inver seRootsToCoef ficients} 

To form the coefficients of a polynomial from its roots is a simple matter. To 
unravel the roots from the coefficients is generally far more difficult. The topic of 
polynomial factorisation is pursued is a subsequent chapter where practical methods 
are presented for extracting the roots of polynomials of high degrees. 

Matrices and Polynomial Algebra 

So far, in representing time-series models, we have used rational polynomial 
operators. The expansion of a rational operator gives rise to an indefinite power 
series. However, when it comes to the numerical representation of a model, one is 
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constrained to work with finite sequences; and therefore it is necessary to truncate 
the power series. Moreover, whenever the concepts of multivariate statistical anal- 
ysis are applied to the problem of estimating the parameters of time-series model, 
it becomes convenient to think in terms of the algebra of vectors and matrices. For 
these reasons, it is important to elucidate the relationship between the algebra of 
coordinate vector spaces and the algebra of polynomials. 



Lower-Triangular Toeplitz Matrices 

Some of the essential aspects of the algebra of polynomials are reflected in the 
algebra of lower-triangular Toeplitz matrices. 

A Toeplitz matrix A — [a,,] is defined by the condition that o,; ; = ai+k,j+k for 
all i,j and k within the allowable range of the indices. The elements of a Toeplitz 
matrix vary only when the difference of the row and column indices varies; and, 
therefore, the generic element can be written as o^- = o,_j . The n x n matrix 
A = [a.i-j] takes the following form: 







a 0 


Oi-l 


0-2 • 


. . Oi_. 






d i 


ao 


0-1 . 


• • 02- 


(2.80) 


A = 


d-2 


d i 


«0 • 


• ■ O3- 






Q-n— 1 


CXn- 2 


On— 3 * 


. . Oq 



A lower-triangular Toeplitz A = [cti--j] has oti-j = 0 whenever i < j. Such a 
matrix is completely specified by its leading vector a = {oo, . . . , o n -i}. This vector 
is provided by the equation a = Ae o where eo is the leading vector of the identity 
matrix of order n which has a unit as its leading element and zeros elsewhere. 
Occasionally, when it is necessary to indicate that A is completely specified by a, 
we shall write A = A (a). 

Any lower-triangular Toeplitz matrix A of order n can be expressed as a 
linear combination of a set of basis matrices /, L, . . . , L" -1 , where the matrix 
L = [ei, . . . , e n _ 2 , 0], which has units on the first subdiagonal and zeros elsewhere, 
is formed from the identity matrix I = [eo, ei, . . . , e n _i] by deleting the leading 
vector and appending a zero vector to the end of the array: 



(2.81) 



L = 



' 0 


0 


0 


... 0 


0 ' 


1 


0 


0 


... 0 


0 


• O 


1 


0 


... 0 


0 


• ‘ O 


0 


0 


... 0 


0 


Lo 


0 


0 


... 1 


0 _ 



This is a matrix analogue of the lag operator. When q < n, the matrix L q , which 
is the gtlr power of L, has units on the gth subdiagonal and zeros elsewhere. When 
q > n the matrix L q is null; and therefore L is said to be nilpotent of degree n. 
Thus the lower-triangular Toeplitz matrix A may be expressed as 



(2.82) 



A — Oo-f a\L -(-•••-)- ot n —iL n 
= Oi(L). 
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This can be construed as polynomial whose argument is the matrix L. The notation 
is confusable with that of a polynomial in the lag operator L which operates on 
the set of infinite sequences. Distinctions can be made by indicating the order of 
the matrix via a subscript. The matrix L a Q is synonymous with the ordinary lag 
operator. 

According to the algebra of polynomials, the product of the ptli degree polyno- 
mial a(z) and the fcth degree polynomial j3{z ) is a polynomial 7(2) = a(z)/3(z) = 
f3(z)a(z) of degree p+k. However, in forming the matrix product AB = a(L)/3(L) 
according the rules of polynomial algebra, it must be recognised that L q = 0 for all 
q > n; which means that the product corresponds to a polynomial of degree n — 1 
at most. The matter is summarised as follows: 

(2.83) If A = a(L) and B = /3(L) are lower-triangular Toeplitz matrices, then 

their product T = AB = BA is also a lower-triangular Toeplitz matrix. 
If the order of T exceeds the degree of 7 (z) = a(z)(3(z) = /3(z)a(z), 
then the leading vector 7 = Tei contains the complete sequence of the 
coefficients of 7 (z). Otherwise it contains a truncated version of the 
sequence. 



If the matrices A, B and T were of infinite order, then the products of multiplying 
polynomials of any degree could accommodated. 

The notable feature of this result is that lower-triangular Toeplitz matrices 
commute in multiplication; and this corresponds to the commutativity of polyno- 
mials in multiplication. 

Example 2.3. Consider the polynomial product 

(2 84) = ( a ° + aiZ + «2* 2 )(A) + Piz) 

= CloA) + (<*oAl + a lPo) z + (ouAl + d 2 A> )z 2 + O^AlA 3 ■ 



This may be compared with the following commutative matrix multiplication: 



(2.85) 



where 



a 0 


0 


0 


0 


Oil 


do 


0 


0 


O' 2 


Oil 


«0 


0 


0 


Ol2 


OL l 


do 



A) 0 0 0 
Ai A) 0 0 
0 Pi A> 0 

0 Oftft 



7o 


0 


0 


0 


7i 


7o 


0 


0 


72 


7i 


7o 


0 


73 


72 


7i 


7o 



(2.86) 



70 = a 0 f3 0 , 

71 = a oPi + ouA>> 

72 = ouAl +a 2 A), 

73 = o^Al- 



The inverse of a lower-triangular Toeplitz matrix A = a(L) is defined by the 
identity 



(2.87) 



A~ 1 A = a~ 1 (L)a(L) = I. 
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Let a~ 4 {z) = {w 0 + 0 J 1 Z+ ■ ■ ■ +co n -iZ n ^ 1 + • • •} denote the expansion of the inverse 
polynomial. Then, when L in put in place of z and when it is recognised that 
L q = 0 for q> n, it will be found that 

(2.88) A — T to\L -(-*** T uj n —iL 
The result may be summarised as follows: 

(2.89) Let a(z ) = ao + aiz + ■ ■ ■ + a p z p be a polynomial of degree p and let 

A — a(L) be a lower-triangular Toeplitz matrix of order n. Then the 
leading vector of A -1 contains the leading coefficients of the expansion 
of a~ 4 (z) = {w 0 + viz + f ^>n-iz n ~ 1 + ••■}. 

Notice that there is no requirement that n > p. When n < p, the elements of the 
inverse matrix A are still provided by the leading coefficients of the expansion 
of a~ 1 (z), despite the fact that the original matrix A = a(L) contains only the 
leading coefficients of a(z). 

Example 2.4. The matrix analogue of the product of 1 — 02 and (1 — Oz ) _1 = 
{1 + 9z+ 9 2 z 2 + ■ ■ •} is 



1 


0 


0 


0 ' 




1 


0 


0 


0 ' 




' 1 


0 


0 


0 ' 


-9 


1 


0 


0 




9 


1 


0 


0 




0 


1 


0 


0 


0 


- 9 


1 


0 




e 2 


9 


1 


0 




0 


0 


1 


0 


0 


0 


-9 


1 




9 3 


9 2 


9 


1 




0 


0 


0 


1 



This matrix equation is also the analogue of the product of (1 — 6z)/( 1 — 9 4 z 4 ) and 
l + 02 + 0 2 2 2 + 6> 3 2 3 . 



Circulant Matrices 

A circulant matrix is a Toeplitz matrix which has the general form of 



(2.91) A = 



a 0 


^n— 1 


a n-2 


. . . Qq 


Oil 


ao 


1 


. . . Ct 2 


012 


a 1 


ao 


. . . «3 


Oin— 1 


Otn—2 


®-n— 3 


. . . ao _ 



The vectors of such a matrix are generated by applying a succession of cyclic 
permutations to the leading vector, which therefore serves to specify the matrix 
completely. The elements of the circulant matrix A = [a,,,'] fulfil the condition that 
a.ij = a{(i — j) mod n}. Hence, the index for the supradiagonal elements, for which 
1 — n < i — j <0, becomes (i — j ) mod n = n + (i — j). 

Any circulant matrix of order n can be expressed as a linear combination of a 
set of basis matrices I,K, . . . , K n ~ 4 , where K = [ei, . . . , e„_i, eo] is formed from 
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the identity matrix / = [eo, e\, . . . , e„_i] by moving the leading vector to the back 
of the array: 



(2.92) 



' 0 


0 


0 . 


. 0 


1 ' 


1 


0 


0 . 


. 0 


0 


• 0 


1 


• 0 


• 0 


• 0 


0 


• • 0 


0 . 


. 0 


0 


0 


0 


0 . 


. 1 


0 



This is a matrix operator which effects a cyclic permutation of the elements of any 
(column) vector which it premultiplies. Thus, an arbitrary circulant matrix A of 
order n can be expressed as 



(2.93) 



A — ocq! T a\K + • • • + cx n — i A n 
= a(K). 



The powers of K form an n-periodic sequence such that K 3+n = K 3 for all 
j or, equivalently, K } j = K ) (j mod n) . The inverse powers of the operator 
K are defined by the condition that K~ q K q = K° = I. It can be confirmed 
directly that K~ q = K n ~ q . However, this also follows formally from the condition 
that K n = K° = I. It may also be confirmed directly that the transpose of K is 
K' = K n ~ 1 = A'" 1 . 

It is easy to see that circulant matrices commute in multiplication, since this 
is a natural consequence of identifying them with polynomials. Thus 



(2.94) If A = a(K) and B = P(K) are circulant matrices, then their product 

T = AB = BA is also a circulant matrix whose leading vector 7 = 
T eo contains the coefficients of the circular convolution of the leading 
vectors a = Ae 0 and (3 = Be q. 



Example 2.5. Consider the following product of circulant matrices: 



(2.95) 



Here 



(2.96) 



CKO 0 ^2 




Po 0 @2 Pi 




7o 73 72 7i 


CKi CKo 0 CK 2 




Pi Po 0 P2 




7i 7o 73 72 


CK 2 CK 1 CKo 0 




P2 Pi Po 0 




72 7i 7o 73 


0 CK 2 CKi CKq 




0 P2 Pi Po _ 




.73 72 7i 7o _ 



70 = a oA) + 0,2^2, 

71 = a i/3o + oto(3\, 

72 = cc2 /3 o + ot \( 3 \ + ao/?2, 

73 = «2/3i +ai/3 2 , 



represent the coefficients of the circular convolution of {ao; a i> a 2, 0 } & n d 
{/3o,f3i,/32,0j. Notice that, with /3 2 = 0, the coefficients { 70 , 71 , 72 , 73 } would 
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be the same as those from the linear convolution depicted under (2.85). Thus 
it is confirmed that the coefficients of the linear convolution of { 0 : 0 , 01 , 0 : 2 } and 
{/3o,/?i} may be obtained by applying the process of circular convolution to the 
padded sequences { 00 , 01 , 02 , 0 } and {/?o, /?i, 0, 0}. 

If A = a (K) is a circulant matrix, then its inverse is also a circulant matrix 
which is defined by the condition 

(2.97) A- 1 A = a~\K)a(K) = I. 



If the roots of o (z) = 0 lie outside the unit circle, then coefficients of the 
expansion a (z)^ 1 = {loq + wiz + ■ ■ ■ + uj n -iz n ^ 1 + ■ ■ •} form a convergent sequence. 
Therefore, by putting K in place of z and noting that K ] q = K } (q mod n), it is 
found that 



(2.98) 



j=0 



3 = 0 



A 1 - + 1 51 W (j"+!) \ K H j ^2 ^(jn+n-1) \K n 1 



1=0 



— 0o + ^iK + • • • + ij) n —iK n 



Given that u>j — > 0 as j — > 00 , it follows that the sequence {"00, ipi, ■ ■ ■ > ipn-i} 
converges to the sequence {wo, wi, . . . , Lc n -i} as n increases. If the roots of 
a(z) = 0 lie inside the unit circle, then it becomes appropriate to express A as 
A = K~ 1 (a n - 1 + a n - 2 K~ x + • • • + ai/\ 2 -" + a 0 I< l ~ n ) = K^a'iR- 1 ) and to 
defined the inverse of A by the condition 

(2.99) A- 1 A = c/- 1 (K- 1 )a , (K- 1 ) = I. 

The expression under (2.98) must then be replaced by a similar expression in terms 
of a convergent sequence of coefficients from the expansion of o:\z~ 1 ). 



The Factorisation of Circulant Matrices 

The matrix operator K has a spectral factorisation which is particularly useful 
in analysing the properties of the discrete Fourier transform. To demonstrate this 
factorisation, we must first define the so-called Fourier matrix. This is a symmetric 
matrix U = n~^ 2 [W jt - 1, j = 0, . . . , n — 1] whose generic element in the jth row 
and tth column is W A = exp(—i2ntj/n). On taking account of the n-periodicity 
of W q = exp(— i2 / nq/n), the matrix can be written explicitly as 

11 1 ... 1 

1 W W 2 ... W 71 " 1 

1 w 2 w 4 . . . w n ~ 2 

1 w n ~ 1 w n ~ 2 ... w 

The second row and the second column of this matrix contain the nth roots of 
unity. The conjugate matrix is defined as U = n~ x ! 2 \W~A^ t, j = 0 , . . . ,n — 1]; 



(2.100) 



U= -j= 
vn 
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and, by using W q = W n q , this can be written explicitly as 

11 1 ... 1 
1 W n ~ 4 W n ~ 2 ... W 

1 w n ~ 2 w n ~ 4 ... w 2 
i w w 2 ... r - 1 

It is readily confirmed that U is a unitary matrix fulfilling the condition 
(2.102) UU =UU = I. 



(2.101) 



U=^= 

\Jn 



To demonstrate this result, consider the generic element in the rtlr row and the stir 
column of the matrix UU = [5 rg ]. This is given by 



(2.103) 



Here there is 






1 

n 



n— 1 

y w rt w~ st 

t - 0 



1 

n 



n— 1 

w (r ~ s)t . 

t - 0 



(2.104) 



= W qt = exp(— i2Trqt/n) 

— cos(— i2Tvqt/n) — ism(—i2Trqt/n). 



Unless <7 = 0, the sums of these trigonometrical functions over an integral number 
of cycles are zero, and therefore W qt = 0. If q = 0, then the sine and cosine 
functions assume the values of zero and unity respectively, and therefore ^ 2 t W qt = 
n. It follows that 



( 1, if r = .s: 

to, if r ± s, 



which proves the result. 

Example 2.6. Consider the matrix 



(2.105) 





'1 1 1 1 




'1 1 1 1 


1 


1 w w 2 w 3 


1 


1 w w 2 w 3 


2 


1 w 2 w 4 w 6 


" 2 


1 W 2 1 w 2 




1 w 3 w 6 w 9 




1 w 3 w 2 w 



The equality comes from the 4-period periodicity of W q = exp(— nq/2). The con- 
jugate matrix is 



(2.106) 





'1 1 1 1 




'1 1 1 1 


1 


1 w ~ 4 w ~ 2 w ~ 3 


1 


1 w 3 w 2 w 


2 


1 w ~ 2 w ~ 4 w ~ 6 


- 2 


1 W 2 1 w 2 




1 w ~ 3 w ~ 6 w ~ 9 




1 w w 2 w 3 
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With W q = exp(— nq/2) = cos(— nq/2) — i sin(— nq/2), it is found that W° = 1, 
W 1 = —i, W 2 = —1 and W 3 = i. Therefore, 



(2.107) UU 





'1 1 1 1 ' 




'1 1 1 1 ' 




' 1 


0 


0 


0 ' 


1 


1 — i — 1 i 




1 i — 1 — i 




0 


1 


0 


0 


4 


1-11-1 




1-11-1 




0 


0 


1 


0 




1 i — 1 — i 




1 — i — 1 i 




0 


0 


0 


1 



Consider postmultiplying the unitary matrix U of (2.100) by a diagonal matrix 







'1 


0 


0 


...o' 






o r- 1 


0 


... 0 


(2.108) 


D = 


0 


0 


w n ~ 2 


... 0 






0 


0 


0 


... w 



Then it is easy to see that 
(2.109) 



UD = KU , 



where K is the circulant operator from (2.92). From this it follows that K = UDU 
and, more generally, that 



(2.110) K q = UD q U. 

By similar means, it can be shown that K' = UDU, where 

(2.111) D = diag{l, W, W 2 , . . . , W" -1 } 



is the conjugate of D. The following conclusions can be reached in a straightforward 
manner: 



(2.112) If A = a(K) is a circulant matrix then 

(i) A = a{K) = Ua{D)U, 

(ii) A' = a(K') = Ua{D)U , 

(iii) A- 1 = a(K ) = Ua~ 1 (D)U. 

(2.113) If the elements of the circulant matrix A are real numbers, then the 
matrix is its own conjugate and 

A = A = Ua(D)U. 



Notice that the set of matrices D k \ k = 0, ...,n — 1 forms a basis for an 
n-dimensional complex vector space comprising all diagonal matrices of order n. 
Therefore, provided that its coefficients can be complex numbers, the polynomial 
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a(D) in the expressions above stands for an arbitrary diagonal matrix with real or 
complex elements. If the coefficients of a(D) are constrained to be real, then the 
jth element of the diagonal matrix takes the form of 

(2.114) Sj = UjW 3 * = Y a ii co s ( w ^) ~ * sin K'*) } - 

3 3 

where uij = 2Trj/n. In that case, the sequence of complex numbers {5j',j = 
0,1, . . . ,n — 1} consists of a real part which is an even or symmetric function of t 
and an imaginary part which is an odd or anti-symmetric function. 

Example 2.7. Consider the equation 

(2.115) X = Ux{D)U = Ux(D)U , 

which defines the real-valued circulant matrix X. The second equality follows from 
the fact that matrix is its own conjugate. Observe that, if eg = [1,0,..., 0] 7 is the 
leading vector of the identity matrix of order n, then 

Xeo X [*^ 0 ? • j l] j 

Ue 0 = Ue 0 = i = [1,1,..., I] 7 , 

(2.116) 

x{D)i = £ = [£o, £i> ■ ■ • , £n-i] and 
x(D)i = £* = [£n-l)£n-2, ■ ■ ■ , ^’o] / - 

Here the vector £ is the discrete Fourier transform of the vector x. Its elements 
are the values {£k = x{zk)\ k = 0 ,n — 1} which come from setting z = zt = 
exp(— 27 t k/n) in x{z) which is the z-transform of the sequence. {xq,Xi, . . . , x„_i}. 
Premultiplying the equation X = Ux(D)U from (2.115) by U and postmultiplying 
it by eo gives 



(2.117) Ux = £, 

which represents the direct Fourier transform of the vector x. Postmultiplying the 
equation by eg gives 

(2.118) x = U£; 

and this represents the inverse Fourier transform by which x is recovered from £. 



Example 2.8. Consider the multiplication of two circulant matrices 



(2.119) 



A = a(K) = Ua(D)U and 
B = a(K) = UP{D)U. 



Their product is 



(2.120) 



AB = Ua(D)UU f3(D)U 
= Ua(D)/3(D)U. 
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On the LHS, there is a matrix multiplication which has already been interpreted 
in Example 2.5 as the circular convolution of the sequences {ace • • • , a n -i} and 
{/3q, . . . ,p n - 1 } which are the coefficients of the polynomials of a(z) and (3(z). On 
the RHS there is a matrix multiplication a{D)(3{D) which represents the pairwise 
multiplication of the corresponding nonzero elements of the diagonal matrices in 
question. These diagonal elements are the values {a(zk);k = 0 , ...,n— 1} and 
{/3(zk)\ k = 0, . . . , n — 1} of the discrete Fourier transforms of the sequences; and 
they come from setting z = Zk = exp(— 2irk/n) in a(z) and j3{z). Thus it can 
be demonstrated that a convolution product in the time domain is equivalent to a 
modulation product in the frequency domain. 
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CHAPTER 3 



Rational Functions 
and Complex Analysis 



The results in the algebra of polynomials which were presented in the previous 
chapter are not, on their own, sufficient for the analysis of time-series models. 
Certain results regarding rational functions of a complex variable are also amongst 
the basic requirements. 

Rational functions may be expanded as Taylor series or, more generally, as 
Laurent series; and the conditions under which such series converge are a matter 
for complex analysis. 

The first part of this chapter provides a reasonably complete treatment of the 
basic algebra of rational functions. Most of the results which are needed in time- 
series analysis are accessible without reference to a more sophisticated theory of 
complex analysis of the sort which pursues general results applicable to unspecified 
functions of the complex variable. However, some of the classic works in time-series 
analysis and signal processing do make extensive use of complex analysis; and they 
are liable to prove inaccessible unless one has studied the rudiments of the subject. 

Our recourse is to present a section of complex analysis which is largely self- 
contained and which might be regarded as surplus to the basic requirements of time- 
series analysis. Nevertheless, it may contribute towards a deeper understanding of 
the mathematical foundations of the subject. 

Rational Functions 

In many respects, the algebra of polynomials is merely an elaboration of the 
algebra of numbers, or of arithmetic in other words. In dividing one number by 
another lesser number, there can be two outcomes. If the quotient is restricted to 
be an integer, then there is liable to be a remainder which is also an integer. If no 
such restriction is imposed, then the quotient may take the form of a interminable 
decimal. Likewise, with polynomials, there is a process of synthetic division which 
generates a remainder, and there is a process of rational expansion which is liable 
to generate an infinite power-series. 

It is often helpful in polynomial algebra, as much as in arithmetic, to be aware 
of the presence of factors which are common to the numerator and the denominator. 

Euclid’s Algorithm 

Euclid’s method, which is familiar from arithmetic, can be used to discover 
whether two polynomials a(z) = a^+ai z+- ■ ■ +a p z p and /3(z) = (3o+f3iZ+ ■ ■ ■ +(3kZ k 
possess a polynomial factor in common. Assume that the degree of a(z) is no less 
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